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ABSTRACT 


ine immediate attention of the control systems engineer 
is directed to the dynamic behavior of the system under 
study. It is important to study the effects on overall 
system performance of varying one or more parameters (mass, 
inertia, gain, resistance, etc.). It is equally important 
to determine whether a desired dynamic behavior can be 
achieved with any set of values for the parameters--if not, 
redesign is indicated. 

In this thesis a control systems analysis package is 
developed using parameter plane methods. It is an inter- 
active, user-friendly computer aid. Given a characteristic 
equation containing two variable parameters, the output of 
the analysis may be either tabular or graphical, with plots 
of any of the following types: 

1. Constant damping curves as a function of frequency, 
2. Constant frequency curves as a function of damping, 


3. Constant sigma lines (real root lines), 


4. Constant zeta-omega (damping-frequency) curves. 


TABLE OF CONTENTS 


Is INTRODUCTION ------- =s SSE 2 
ІІ. DERIVATION OF PARAMETER PLANE EQUATIONS-- -- © 
III. APPLICATIONS OF THE PARAMETER PLANE METHOD------- 20 
A. ALGEBRAIC SOLUTION--------------------------- 20 

1. Feedback Compensation-------------------- 20 

2. Cascade Compensation--------------------- Be 

В. DOMINANCY OF THE SPECIFIED ROOTS------------- 38! 

IV. PROGRAM DESCRIPTION------------------------------ 50 
А. ТНЕ РЮКОСВАМ---------------------------------- 50 

1. Constant Zeta Contours------------------- ou 

2. Constant Omega Contours------------------ ӨӨ 

3. Constant Sigma Contours------------------ 58 

4. Constant Zeta-Omega Curves--------------- 54 

о. Tabular Output--------------------------- во 

6. Plot Кезса1іпсе--------------------------- 56 

B. INSTRUCTION TO THE USER---------------------- one 

V. PARAMETER PLANE CURVES-GRAPHICAL METHOD---------- 61 
A. GRAPHICAL SOLUTION--------------------------- 61 


B. MISCELLANEOUS ASPECTS OF THE PARAMETER PLANE- 87 


VI. CONCLUSIONS =- === Со ТЕБЕ: SSE PE s 90 
APPENDIX A - FUNCTIONS U, (§)--------------------------- 92 
APPENDIX B - PARAMETER PLANE PROGRAM SOURCE CODE------- 93 
Ι.15τ ΟΕ ΒΕΕΕΒΕΝΕΕ5-.---------..----....--.-- ел ишш 11: 
INITIAL DISTRIBUTION LIST- -a a --"Ἠπὶ 115 


I. INTRODUCTION 


The analysis and synthesis of linear feedback control 
systems, or the compensation of same, can be realized by 
three general methods. The first of these can be called the 
integral method. Given a control system, described by a set 
of differential equations, one selects a cost function to be 
minimized with respect to certain variable system parameters. 
The major drawback with this method is the difficulty of 
varying more than one parameter at a time. The second 
method is the Bode frequeney response technique whereby the 
system's open loop transfer function is manipulated to obtain 
the desired system response. This method also has its 
inherent weaknesses: difficulty of application to non-unity 
feedback control systems, difficulty in interpreting the 
closed loop transient response in terms of the open loop 
frequency response, and difficulty of varying more than one 
parameter. Third are the algebraic methods. Within this 
category can be included the familiar root locus method. 
Here, a graphical technique is provided by which the set of 
all points which could potentially be made roots are plotted 
in the S-plane. The root locus method is a valuable and 
powerful tool when only one parameter is varied; results 
are less satisfactory for two parameters and of little use 
when three or more parameters are involved. 
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Methods for studying the parameter-root relationship 
when two or more parameters are variable are clearly of 
considerable value. For a linear system, the set of 
differential equations that describe that system can be 
transformed into algebraic equations and manipulated to 
provide a characteristic polynomial. Since the coefficients 
of the characteristic polynomial are deterimined by the 
system parameters, ác follows that some relationship exists 
between the value of any parameter and the value of the 
characteristic roots. In reference (1), Mitrovic developed 
an algebraic/graphical method for obtaining the roots of a 
polynomial in terms of two variable parameters. In 
references (2), (3), and (4), Choe, Hyon, and Nutting, 
respectively developed and extended the Mitrovic method to 
the compensation of linear continuous feedback control 
systems. The disadvantage of the Mitrovic method is that the 
variable parameters may appear in no more than two coeffi- 
cients of the characteristic equation, which limits the 
flexibility of the technique. In reference (5), Siljak 
introduced a method for obtaining the roots of a polynomial 
in terms of two variable parameters that may appear in any 
and all of the coefficients of the polynomial. Later, x 
and Towill [Ref. 6] extended this method to the compensation 
of linear continuous feedback control systems. It is from 


the latter work that the ensuing parameter plane equations 


were developed. General methods of compensation will be 
presented, and an attempt will be made to relate the root 
locus and the parameter plane methods as a set of comple- 
mentary techniques which, when applied in tandem, represent 
the most satisfactory tool to date for designing linear 
feedback control systems. 

The parameter plane method, which works well for two 
variable parameters and which may be extended to three or 
more parameters, is purely algebraic, and the resulting plots 
are valuable aids to analysis. The term parameter plane 
comes from the plot for two parameters--in a rectangular 
coordinate space one parameter will define the abscissa while 
the second parameter defines the ordinate (the S-plane is 
inconvenient for presenting the desired results). Three 
parameters define a 3-dimensional parameter space, etc. For 
design problems it is convenient to think of the algebraic 
calculations as a mapping procedure. By choosing a point on 
the S-plane, the characteristic polynomial acts as a mapping 
function whereby the point may be "mapped" onto the alpha- 
beta plane (alpha and beta are the two variable parameters 
to be used throughout the remainder of this text). The 
relationship between being able to place the roots of a 
polynomial at specific locations in the S-plane and the 
compensation of linear feedback control systems is as 


follows. А feedback control system, including any added 


compensators which may contain variables, can be reduced to a 
ratio of two polynomials (the closed loop transfer function). 
A specified system response in terms of overshoot, bandwidth, 
settling time, steady-state accuracy, etc., can theoretically 
be obtained by placing a pair of complex conjugate roots of 
the characteristic equation at a specific location in the S- 
plane, while ensuring that the real part of this complex root 
pair (the dominant roots) is smaller in magnitude than the 
real parts of the remaining roots of the characteristic 
equation. The problem of compensation, thus of feedback 
control system design, reduces to one of placing the dominant 
roots of the characteristic equation at the desired location. 
The ability of the parameter plane method to achieve this 


goal will become obvious. 


ІІ. DERIVATION OF PARAMETER PLANE EQUATIONS 


A linear feedback control system's characteristic equation 


can be expressed as a polynomial of the following form: 





E k 
#(5) + 5 aV S - 0, where (2-1) 
k=0 
a. (k20,1,...,m) are real coefficients 
S = -otjw = -Ewt jw Dre 


ω is the undamped natural frequency and 
50 is the relative damping coefficient 
In reference (5) it is noted that 58 ша ὃς тергезепгеа Бу 


the following: 


Eo 0D (- 53v 1-£^ 9 (- 6) (2-2) 


where 


М k " k+l. ,, 
T(-£)=(-1) "7, () and U. (-£)=(-1) 710 (8). 


T, C&) and U, CS) are Chebishev functions of the first and 
second kind respectively. Values of zeta and omega will be 
Eoussdered such that O«£«1 and O«o«e. Values of Т, апа Ше 

are tabulated in various appendixes. More important to digital 


computer analysis, they can be obtained from the following 


recursive relations: 


lI 
O 


DC SEDED се 


(2-3) 


| 
© 


КИС ee 


Here, Tg(C5)71, Т/(%)=Ёё, 0258250, г, (52-1. SubstTtutaws 
equation (2-2) into (2-1) and setting the real and imaginary 


parts to zero independently, one obtains: 


hws 


k = 
a, 9 T, (-£)-0 


k=0 


(2-4) 
т 
k 
> а Т.(-8)-0 
k=0 k k 


Employing equations (2-3), one obtaihs from equation (2-4): 


(-1) а ши (Е) - ο 


има 


k=0 


(2 


AC i 
Е (-1) а 0 U, C5) - 0 


IMS 


k 


Now consider the coefficients а, of the characteristic 
equation (2-1) as linear functions of the variable system 


parameters, o and f, as follows: 


а = bw + св * d, (2-6) 
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C ils this relation for a equations (2-5) become: 


ка 


а В. + BC, + Di Z 
E2 a) 
! B = 
ob, + Со + D, 0 
where 
m m 
Ву = 5 (-106 вт у B,- z (-D* bou 
к=) К=0 
m m ў 
C, - > (-1)* cet u у Со = z (-D* ew Uu 
110, [<= 0 
(2-8) 
m m 
MEX (Daetu. DNE M MENU 
k=0 = k=0 


Since equations (2-7) are linear in the two unknowns’ alpha and 


beta, Cramer's rule may be applied to obtain: 


ΝΕ απ» 


TUE пре лири. (2-9) 
B1C5-B4C, В,62-В264 


л. 


Equations (2-9) are now functions of zeta and omega. Hence, 
by fixing either zeta or omega and varying the remaining 
parameter, the constant omega or constant zeta S-plane 
contours respectively can be mapped into the real domain 
of the alpha-beta or parameter plane. 


In reference (5) the following relationships are noted: 





9 
= 4 1 у = 
S D jov 1-£ Q. 


2 
k+1 k k-1 


оу + 0 (2-10) 


rÜ 
-= 
| 
Ц 
t 
Str 
© 
- 
€ 
N 
&) 
T 
j= 


Р, апа Qi are related to the Chebishev functions by: 


“wk M мек 
ED T, C-£) = (-1) & T. š) 


(2-11) 


< иЕЦу (-Е) # (ро 15 


< 
= 
I 


(2) 
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By employing equations (2-10) and (2-11), one obtains 


(proceeding as before); 


m m 
> аю = 0 XO - 0 ο ιδ) 


Combining equations (2-6) and (2-12) with Cramer's rule, one 


again arrives at equations (2-9) where the following expressions 


now apply: 
m m 
Е > D. Q = 0 отук = 0 
1 0 К "k-1 2 ET κ “k 
m m 
E —» c, Q = 0 = oc, OG = 0 (2-13) 
1 E κ “k-1 2 k-0 κ “k 
m m 
m= 2 d, Q = 0 5 ο 
1 к-0 k “k-1 2 k-0 κ Kk 


Equations (2-9) and (2-13) are useful for mapping constant 
zeta-omega curves from the S-plane into the parameter plane. 
As will be demonstrated later, these curves play an important 
role in dominance considerations. 

If the complex variable S is substituted in equation (2-1) 


by letting S -—o, where sigma corresponds to values of S along 
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the real axis, then according to equation (2-6) the 


characteristic equation (2-1) becomes: 


m 
(-1)*b, 0% + Bs (=l) e g 


m 
Q D 
=0 k=0 k=0 


k 


(2-14) 


The above expression represents a straight line in the alpha- 
beta plane for a given value of sigma. Hence a point on the 
real axis in the S-plane maps into a straight line in the 
alpha-beta plane. In addition, for given values of alpha, 
beta, and sigma which satisfy equation (2-14), the 
characteristic equation (2-1) must have a real root at minus 
sigma. On the constant zeta and omega curves previously 


defined, for certain values of alpha and beta (say, for 


values obtained from equations (2-9) for given values of 


zeta and omega) the characteristic equation will have a pair 





от сопріех σοπ]μσαυς Ἱ9θουδθ r = eG ju mw В 

The significance of the above discussion is that by 
applying equations (2-9) and (2-14) one can, for a specified 
value of zeta, omega, and Sigma, compute the value of alpha 


and beta such that the characteristic equation will have a 
2 
1 1 i 


remaining roots of the characteristic equation can then be 


pair of roots at S - -Е a jo 1-£ το nl 


calculated by dividing out the two known or specified 


roots. This method, where zeta, omega and sigma, or simply 
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zeta and omega are specified, and where the computations for 
аа апа beta are done algebraically, will be referred to as 
ШИЕ algebraic parameter plane solution. 

To solve the problem in general for all values of zeta, 
omega, and sigma, it becomes necessary to plot a family of 
parameter plane curves for various values of zeta, omega, 
sigma, and if desired, zeta-omega. On the resulting 
parameter plane plot one can, by choosing an operating point, 
graphically read from the curves the values of alpha and beta 
and, hence, the values of the m roots of the corresponding 
"m order characteristic equation. This latter method will 
be referred to as the graphical parameter plane solution. 

The algebraic solution has the advantage that the labor 
of plotting the curves сап be avoided, but the disadvantage 
remains that without the curves it is difficult to pick the 
Optimum values of zeta and omega so as to ensure dominance 
while still meeting the system specifications. The graphical 
solution has the advantage that one has a "picture" of the 
way the characteristic roots move about in the S-plane as 
alpha and beta are varied. This enables one to choose the 
values of alpha and beta corresponding to the best values of 
zeta, omega, sigma, and zeta-omega for all roots of the 
E uccteristic equation. This feature of the parameter 
пе points out a strong justification for attempting to 


obtain the parameter plane curves. And with the employment 


15 


of a digital computer and an appropriate algorithm to realize 
the parameter plane curves, the advantage of the algebraic 
method becomes muted. 

Recursion methods (equation (2-3)) are by no means the 
only methods of producing algebraic and graphical parameter 
plane data. Thaler and Karmarkar [Ref. 7] describe a matrix 
solution to the parameter plane problem. Essentially, a 
matrix of coefficients may be manipulated to obtain the 


following general form: 


Cy -W 0 0 0 di e а 
т а ы т 0 а е 8 
2 vedi 2 2 

m-2 

-1 -2£o В. О В 

жи т-2 
-1 ала - : : l -0 
m-2 
2 
= 41-2 2 
2 

-2#о d 217% d 

пе = а 
ca О 0 0 1 d, 2ξω е В 
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where б,» сұ» d, w, and w are as described before, and: 


ШЕШ m) аге the coefficients associated with the 


e 
E non-linear alpha-beta product terms 


€ 


ПІ--< 
RH. is the sum of the m-2 roots of the polynomial character- 
istic equation taken one at a time 


В. 1-9 ο πο оет =лтосошы taken m-2 at a time 


u cr sby the application of appropriate row operations, 


this matrix may be reduced to the following row-echelon form: 





О О О Куу Kio α 

1 0 О Kod Коо В 
m-2 

О 1 О Και Koo Ж, 
= Q 
oes) 
1 
ав 
О О О 1 k nl E 





T 7 


For the case when all coefficients of the characteristic 
equation are linear, i.e. ey, Ck71,...,m) = 0, then Κο 
(Ке1,...,п) - ОР "πα 


о- -Қ 


Jl 
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| 
| 
B= -K 
as obtained from the first two rows of the matrix equat Toms x 
One should note that in arriving at equations (2-15), x 
approximately я row operations аге required for the row- 

echelon matrix formulation. for.each point ofthe parameiem x 
plane curves (e.g., each time either zeta or sigma are | 
varied). Compare this with the approximate m calculations 
required to obtain the recursion equations of the previous 

chapter, and the matrix method becomes relatively inefficient 

for larger order systems. 

One should not, however, discard the matrix approach | 
entirely. For small order characteristic equations, thie 
technique compares favorably with the recursion method. 

And when the variable parameters are non-linear--when one 
must deal with alpha-beta product terms--the matrix approach 
affords a more direct method of obtaining the alpha-beta 


2. Whether the recursion or matrix method is utilized 


1 
From equations (2-15), one obtains the two quadratic forms 


2 z = + 
ee (Ко Ко КК. Кро 16 Е 


ШЕ. 
K,,8^*(K,, K,,-K,, К.,2%1)8%К, 150 


from which alpha and beta are easily derived. 
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ur depend on the inclusion of alpha-beta product terms; 


ultimately, it is a matter of personal preference. 


15, 


ІІІ. APPLICATION OF THE PARAMETER PLANE METHOD 


A. ALGEBRAIC SOLUTION 

In this section it will be assumed that the system 
performance specifications can be met by placing a pair of 
complex conjugate roots at a specific location (i.e., by 
choosing appropriate values for zeta and omega). If after 
computation of the necessary values of alpha and beta to 
locate the roots as desired it is found that these specified 
roots are not dominant, then either a different value of zeta 
and/or omega must be used (possibly at the sacrifice of some 
performance "ене ры or a different method of compensation 
will have to be attempted. In a later section a method will 
be addressed whereby the dominancy requirement may be. 
achieved. 

1. Feedback Compensation 


For a unity feedback Controlesystemau et 





Se t K . ай 
G = = Τ (ЭЖИ; 


where K is the forward path gain (a variable) and e(S) is a 
polynomial in S representing the poles of the open loop 


transfer function of the uncompensated system. In equation 
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(3-1), L corresponds to the system type--for a type 0 
EE ІГ-0, Тог type l, L-1, etc. The system's error 
coefficient is defined as: 


"e L 
K = lim 5 G. (3-2) 


550 т 
мһеге Goc is the open loop transfer function of the 
compensated system. 
| а. Tachometer Plus Acceleration Feedback 
In order to achieve the system performance 
Specifications, a feedback compensator must be introduced. 


рет 
nu sk s^ 
t a 


The resulting compensated system's characteristic equation 


becomes: 


24 . 
ЩО) + К(К 5+К 5) + 0 (3-3) 


and by expanding e(S), equation (3-3) becomes: 


m gm-1 


2 " 
— +...+(е„+КК 25 +(e, +KK, )Ste +K = Ü ( 3-4) 


where L is zero for a type O system (the most general case). 
The following results also apply to a type 1 system if e 15 
РБ 50 zero, and, similarly, for а type 2 system if both 

EG and е) аге set to zero, etc. Combining equations (3-2), 


(3-3), and (3-4) the error coefficient becomes: 
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О 
Е 2-2 (3-5) 


е 550 e(S)*K(K, S«K. S^) са 


or for a type l uncompensated system: 


K 


K = ———— (3-6) 
е e, TKK, 


or for a type 2 uncompensated system: 


^X 
cs KK, | Eo 


Note that if the uncompensated system is type 2, the 
compensated system would be type 1 if tachometer feedback 
or tachometer plus acceleration feedback is used. 

In the compensated system's characteristic 
equation (3-4) let alpha = KK and beta = KK,. Equation 


t 
(3-4) then becomes: 


cu + е Smal 


2 =: 
m-1 +...+(е„+а)5 ые + шеке тк = 0 


Recalling equation (2-6) where in general the coefficients 


of the characteristic equation are of the form: 


а, - Бас Ва 


and letting m-k, then from equations (2-8) one obtains: 
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2 7 о - 
` os -^ zu US 
a ou, O os oe (3-8) 
m m 
D, = z (-D'qj"U , Dy = z (1a, 
k=0 k=0 
since 14-0 апа 01-1. From equations (2-9) one derives: 
E k 
БИО ОС р ζ΄ m 
pa» k=0 = Са к=» 
ЕНЕСІ | з |. 2, οσα) αν Ux d 
2 -ω k=0 
(3-9) 
m 
κα. КШ 
В = S а» = 
ν. ΠΤΙ ΜΕ’ 


If alpha and beta are bees functions of6É Ве forward 
path gain, one can use the steady-state error specification 
to define K in terms of alpha and/or beta. Since zeta and 
omega were assumed to be specified, then from equations 
(3-9) one may solve for alpha and beta. From this, K. and 
K. are readily determined. 
Example 3-1 
The system of Figure (3-1) is to be compensated by using 


tachometer plus acceleration feedback. The system 


Specifications are as follows: 
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K 
ο Ως ο) 


= 






Figure 3-1 


1. Complex roots corresponding to £-0.7 and wz10. 


2. K 29 
е— 


From equation (3-2): 





From which K 212* 6KK, . The compensated characteristic 


equation is 


S"«s^(3*KK )*S(2*KK, )*K T (3-10) 
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Letting alpha = ΚΚ, and beta - KK, equation (3-10) becomes: 


57457(34+0)4S(24+8)+K = о Sg 


From equations (3-8): 


В. = 100 В. = 140 
C1 - 0 с. = -10 
D4 ш -1100-К Di = -1120 


and from equations (3-9): 


_ 10(-1100~K) 140(-1100-К ) +56000 


71000 on -1000 (3-11) 


ШЕШ Пе steady-state accuracy specifications, then, it is 
necessary that K>12+68; let К=12+68. From equations (3-11) 
it is found that 82623, hence K=3750. Therefore, @=48.5, 


and since Q-KR and B-KK, : 


K, E 3750 = 0.0129 
оз _ 
K. = 3750 - 0.1661 
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The compensated system's characteristic equation becomes 


53+51.552+6255+3750=0 “ 


Now zeta=0.7 and w=10 corresponds to 52+145+100=0. 
Dividing equation (3-12) by this quadratic, the remainder 
is St37.5. Since zeta-omega of the desired roots = 7««37.5, 
the complex roots are dominant and the problem is solved. 
b. Tachometer Feedback Only 


Let H = KS. The characteristic equation of the 


compensated system becomes: 


1 


m- 2 = 
S +e 15 +...+e,8 +(e, +a )Ste,+8=0 


2 


Proceeding as in the previous example, one obtains: 


В. = 0 В, αμ 

οι - -ᾱ ο ο 
m m 

роо πετ D, = z (-1)*d,0 Шш 
к-0 к-0 
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For a specified value of zeta and omega, alpha and beta can 
Dp i Qqined from equations (3-13). The error coefficient is 
Len determined directly from equations (3-5), (3-6), or 
(3-7). Thus the error coefficient is fixed for a given 
value of zeta and omega, and if this parameter is to be met, 
the values of zeta and omega may require adjustment. One 
possible approach might be to fix zeta at some value, 

whereby from the given K and equations (3-5), (3-6), or 
(3-7) alpha could be computed. Equations (3-13) could then 
be solved for, first, omega and then beta. The calculations 


Would prove tedious, however. 
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Ехатр1е 3-2 
The same system as used in example (3-1) will be studied 
here, this time with tachometer feedback alone. The same 
system performance specifications are to be met; names 
K 26, zeta = 0.7, and omega = 10. The compensated system's 
characteristic equation becomes: 


59 +357+ (2+КК Б+К т; (3-14) 


Letting а-КК. and B= K here, equation (3-14) becomes: 


53+352(2+0)5+8 = 0 


From equation (3-13) it is found that: 
С - -2%30(1.4)-100(0.96) - -56 


Since alpha is negative, it is seen that positive tachometer 
feedback is required. Further, it is found that the 
remaining root (when equation (3-14) is divided by 
S“+14S+100) is positive; the system is unstable. Hence 
the desired system specifications cannot be met with 
tachometer feedback alone. 

C. Acceleration Feedback Only 

2 


Let H = κ» . The characteristic equation of 


the compensated system becomes: 


б+е +К = 0 
О 


2 
да 4e 
S +e S t... (e, KK. )8 е, 


Proceeding as before, where now а-кк, and B=K: 


= „© . = ω 2 
Б B c а Uo 
SEU. - LS 
m m 
ру = Хх (ш аы - Dona: (-1) а окр 
Кеш k=0 
Solving for alpha and beta yields: 
-D 
2 1 k k-2 
Q = ——_— = - --- == 0) 
D U (-1) d Uk 
Қ Us 2 
OS ST) 
m m 
B= z (-)ба а, 1-2- E (C1f"as"u, 
k=0 2 k=0 р 


Euwsutations for alpha, beta, and K. are performed in the 
same manner as with the preceding tachometer feedback 
example. 

Example 3-3 


The same system of examples (3-1) and (3-2) will now be 


Compensated using acceleration feedback alone. As before, 
К 26, zeta=0.7, and omega-10. Therefore K =2 and the error 
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coefficient is unaffected by the acceleration feedback. Hence 
one can conveniently choose К=12 to meet the ресет 


The compensated system's characteristic equation becomes: 
3 2 
S +(3+KK_)S +2S+K = 0 (3-16) 


If K in equation (3-16) is set equal to 12 as prescribed, only 
one parameter remains and the parameter plane equations 
produce an indeterminate solution. If K is left as the 
variable beta, then equation (3-16) becomes (after the usual 


SubStitutioms): 
s°+(34+a)87+28+B -- 


Ву employing equations (3-15), one obtains @=4 and B= -700. 
since beta is negative, it is concluded that the desired 
roots (i.e., desired values of zeta and omega) cannot be 
realized using acceleration feedback alone, and of course 
neither can the desired error specification be obtained. 
One would therefore choose an alternate method of 
compensation. 

If one chooses to use feedback compensation then 
perhaps tachometer plus acceleration feedback might be 
attempted first using equations (3-9) and the appropriate 
steady-state error specification. If the specifications cannot 


be met in this manner, then it follows that neither 
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En mter nor acceleration feedback alone wi1]1 suffice. ια 
Mee eaS e either the system's specifications must be eased or 
пет type of compensation must be utilized. Pe ites 

ОШО that the specifications are achievable with the combined 
Eu5meter and acceleration feedback, then, if desired, 

ιν ons (3-13) and (3-15) can be employed to investigate 
Eweasibility of tachometer Drei А ИА S, 
ши ресттуету. 


ae Cree Form Whee Feedback Is Not Available Near 
πο ο ναι Path amplifier 


1 





Figure 3-2 


Figure 3-2 shows a system similar to that used 
шинекатр е (3-1) except that now the feedback is inserted at 
E coutput terminals of the amplifier represented by gain К. 
Nu illustrates a system for which it may not be possible 
ENEUrCactical to access the input terminals of the error 


νον. This problem will be solved by means of an example. 
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Ехатр1е 3-4 
As before, the same system specifications are to be met, 
leew. K 26, zeta=0.7, and omega=10. The characteristic 
equation becomes: 
59 (зек )S2+(2+K )S+K = 0 


1 α- = 9 
Letting K and В K. 


59 +(3+0)82+(2+8) 8+К 


II 
c 


(3-17) 


Comparison of equations (3-17) and (3-10a) show that they 
are identical, that is, the solution obtained for alpha and 
beta in example (3-1) applies. There, alpha and beta were 
found to be 48.5 and 623, respectively, while K=3750. For 
the present example no further computations are necessary 


Tomb Lod K and K since they are now the parameters alpha 


t? 
and beta. This points out an important advantage of the 
parameter plane method, namely, that the solutions depend only 
on the characteristic equation and not on the system from 
which the characteristic equation was formed. This principle 
can similarly be applied to control problems involving 
tachometer or acceleration feedback alone. 

2. Cascade Compensation 


For a unity feedback control system let G have the 


form ὋΣ ΘΟΕ ΠΟΠ; 
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К К 
КО = —_— n (3-1) 
EDU c cs" 
where again K is the forward path gain (a variable) and 
e(s) is a polynomial in S representing the poles of the 
τον transfer function of the uncompensated system. 
The letter L again indicates the system type. If in order 


to satisfy the system's requirements a cascade compensator 


ας is required, then let: 


ШЕР STA) 
C Z(S+P) 


G 
Wrth à d.c. gain of unity, placement of this compensator in 
the forward path will not affect steady-state accuracy. 
With Go as indicated here, the values of P and Z are 
computed to obtain the desired system response. If P is 
Tess than 2, a lag network is required and the factor of 
z Of the compensator is inherently present due to the 
physical nature of the compensator (usually an R-C network). 
In this case all forward path amplifier gains can remain 
unaltered to meet the specified accuracy demands. If, 
however, the computed value of P is greater than Z2, a 
lead network is required and the compensated system's 
forward path gain must be raised by the factor of 2 to 
Meer the accuracy specifications. Ав the physical nature 


ОТ the lead network is such that the factor is not 


P 
А 
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inherently present, this factor must be provided eithereby 
adding an amplifier in cascade with the lead network or by 
raising the gain K of the existing amplifier as required to 
achieve steady-state accuracy. 

Continuing, the compensated system's forward path 
Сгап= вети илет олие 
К P S+Z _ y K 


Сес = саб = GU στ ΙΕΡΟΣ “отр е SOEY 


Applying the definition of the error coefficient one 


obtains: 
P 
I ome Y K 
А пори аза 
5-0 L 


and again assuming a type O system where L=0, the compensated 


system's characteristic equation becomes: 


e(s)(S*P)*K Y( S47) ER 


or after expansion: 


m+] m m-1 
+( D+ + + ч 
(Р+е EP (Pe amd NS 


= 1 m-2 


*(Pej*e, )S^«(KY*Pe, *e )S*P(e +K) 5/1 (3-18) 
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Letting ®=P and B=y, equation (3-18) becomes: 


it m m-1 
α α 
E MEDIE ue те 225 iol 


ЕШ” 


2 n 
+(e,ate, )S + (КВ +е а+е „)5+0 (е +K) - (0 (στης) 


| 


Comparison of equation (3-19) with the general form of the 
characteristic equation as specified in equations (2-1) and 


(2-6), it is apparent that K-m-*1, bg-egtk, с 79 -0, b, =e,» 


b =0, а= etc. 


Шы 2 9%. 2 DISC 


It is important to note that the parameter plane 


C =k, d 


Í 
variable beta represents the pole-to-zero ratio of the 
cascade compensator. The S-plane can be divided into regions 
where lag compensation or lead compensation is needed. By 
mapping of variables in the above manner, the parameter 

plane can effectively be divided into corresponding regions 
above and below the straight line 8-1. Then, for values of 
beta less than one a lag network is required and for beta 
greater than one a lead network is needed. In addition, if 
beta is less than 0.1 or greater than 10, a multiple lag or 
multiple lead network, respectively, is required. 

Based on equations (3-19) and (2-8) it is found that: 


2 k-2, k-2 


-(egtK)to саат) ω ου LT 


UJ 
| 


О 


+(-1) - 


"ы © 


Но 


=: --2 ,E-25 . 
D. = o e C 61-2 г, за ӨНЕ S δι ο 
k 
Pree 
(-1) 0,21 
я B2 ко қа k-1 k-1 
В, = "е е202%... .+(-1) a DN -1) 0) Uk. 1 
Со = -—WK 
_ „2 к-2 5-2 
D. = сестра е,02%...%(-1) 9m-2 δι ο 
(3-20) 
k-1 k-lg 
τν ἀν... : 


and from equations (2-9) it is оппа шү: 


_ -D н B,D- 12 
би арс. В - c 


ή απ 


(3-21) 


For a type 1 system, ей in equations (3-20) is set equal to 
Zero, for a type 2 system πα ος etc. On the basis of 
equations (3-20) and (3-21) a cascade compensator can be 
designed. 

Example 3-5 
For the system of Figure (3-3) it is desired to design a 


cascade compensator which places a pair of roots at zeta=0.5 


and omega=1. The error coefficient К should be 50. 
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+ (S+P/ ) 
= S+P 












Figure 3-3 


mers apparent from Figtire (3-3) ttYat K, 10 | 


The characteristic equation is: 
1 3 2 ЕРА 
ENESUSEPOS-*(17*8P)S + (10+17Р+К У)5+Р( 1О+К )<0 
and by the substitutions of 4-P and b= т: 
4 3 2 
ШЕР(8ға)5 ғ(17ғ8о)5 +(1041747+5008)St510¢ = 0 


Applying equations (3-20) and (3-21) one obtains: 


B, - -503 B, = -9 

Е С, - -500 

Ш - 9 D, = 6 

а = 0.0179 = р В = 0.0117 = Y 
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р 
and since Y= 


7» 271.529. This is a lag network for which the 


factor 9.0117 is inherent in the R-C filter design. 

Although treatment will not be presented here, the 
algebraic application of the parameter plane technique can 
be readily applied to combination cascade and feedback 


compensation. 


В. DOMINANCY OF THE SPECIFIED ROOTS 

In the preceding section nothing was done in the 
calculations to make the specified roots a dominant pair. 
As mentioned earlier, the ability to predict a system's 
response on the basis of the location of a pair of complex 
conjugate roots was based on the assumption that the 
magnitude of the real part of the specified roots was much 
less than that of the remaining roots of the characteristic 
equation. In practice, if the real part of the specified or 
primary roots is one half to one fifth or less of the real 
parts of all secondary roots, the system is said to be 
dominant in the primary roots. In many cases the system 
will still meet the specifications even if two pairs of complex 
roots have the same real part, provided the zetas for both 
pairs of roots meet the specifications, and the undamped 
natural frequencies are such that the component time responses 
are not highly additive. Further, even if there exists a 
characteristic root whose real part is closer to the origin 


than that of the primary pair, the presence of a closed-loop 
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zero could make the residue of the close-in root negligible 
as compared to the primary root. If possible, however, one 
ο ον attempts to make the real parts of all secondary 
roots large in magnitude. 

In the preceding examples it should be pointed out that 
in many cases there were actually three and possibly four 
variable parameters. For instance, the forward path gain was 
usually a fixed value in the computations so as to meet the 
minimum steady-state accuracy requirements. There is, 
however, no reason why the gain cannot be raised above the 
minimum value, thus permitting a third degree of freedom. 
When cascade and feedback compensation are used simultaneously, 
the forward path gain and tachometer gain become the third 
and .fourth parameters. 


Recall that the system characteristic equation has the 
m 





ШЕП і(5)- X a = where a,=b,at+c,Bt+d,. To realize the 
| k=0 k κ Kk k k 
System specifications, one places a pair of complex roots 
2 | 2 2 
μα... 10 = + =: 
аф 5--%,691%) y E] >» which implies that S *25.04St*0 7-0 


Since 51 and 9, are known, the quadratic can be divided 


1 
out of the characteristic equation, leaving a polynomial 
which contains the secondary roots of the characteristic 
equation. Since only two of the variable parameters were 
used 1n fixing the primary roots, the remaining parameters 
will appear in the coefficients of the quotient polynomial, 


and it is these parameters that can be varied to achieve 


dominance. 
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Instead of division to obtain the quotient polynomials 
coefficients of like powers will be equated to achieve a. 
system of equations. Let the quotient polynomial be 
given by: 

f(s) = f S” = O (3-22) 
where n=m-2, i.e., equation (3-22) is of order two less than 
the characteristic equation. Applying equations (2-1), 
(3-22), and the quadratic it follows that: 

n 


2 2 ш 
(5 +26, 08+) )¢ X. d S ME es (35-55 
к-0 к= 


Taking а cl and equating coefficients of like power: 


a = f =] 
m n 
а л μι 
а = f E stie ας e 
m-2 п-2 11 oi 1 
(3-24) 

: 2 
аа оча 

2 

= E Е 

ο ο е" 
aoc 
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Equations (3-23) and (3-24) can be solved for the coefficients 
Т іп terms of the coefficients a. The results will be 
applied to the following cases: 

Case of k-3, n-1 


Equation (3-23) becomes: 


s s оға 


2 077 E 
(S^*2&, 0, sca  )( f, Se f9)-8" *a,, {Sta 


Equating coefficients of like power one obtains: 


ae = fi - 1 
а) {02519111 
a. = f 2 2t DET 
1 Шо 
К 2 
ap + 1094 


nin; for the coefficients f results in: 


fi - 1 
2 
80 (а -ω)) 
fo = m = DNE = 2-25-01 (3-25) 
1 1271 
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Case of k-4, n-2 


Proceeding as before, the a coefficients become: 


а, = То - 1 
ag ~  Π 
(3-26) 
2 

= 0) ; 
a M uc MEL. 
a4 = fw, t25,0,f5 
ао = τη 





τ. πι 
е а а зан Е ο А A snn 
€) , 

1 3 1] 1 “©” т 54 2 Q 3 2 1 ü) 

80 а 257 2 2 222 
ο 72 7 Sr.u. 7-9: 1 = аон И 

a) 1 1 

1 
(3-28 
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Case of k=5, n=3 
wmi larly: 


as = Í, = 1 








4 2 
Е с ого ғғ 
D 1 EO 
a τ΄ Έπος 0) 
9 о 255 аан 
= оф уэ; ш г 
ee ii -°1 1-0 
. 2 
8-15 
апа: 
= | 
a 
к: а 2 
f) - a,-2& 5, = 2 
1 
in 73 er 
P = Е 
GQ) 
25,490 2591 
E. 8] 25 180 E 
1 um m2 = 
i 
а, 
ig 7 E 
о) 
1 


въ 





О 





9 
Ў е , 13 и 2 а, 
0524 4 -3 
1 1 1 
а) 1 
4725 
4) 
; 1 
2 2 2 
ο ВЕ оТ 
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Although the coefficientsof f have been derived вог опо 
the fifth order case, they can easily be obtained for haem 


order cases if necessary. 


Example 3-6 (Third Order Characteristic Equation) 


S (S+P/ ) K 
т 5+Р 5(5+4) 





Figure 3-4 


Design a cascade compensator for the system of Figure (3-4) 
to obtain: 

1. Characteristic roots at zeta=0.5 and omega=40. 

2. К 2250. 

e~ 

3. The specified roots are to be dominant. 

The characteristic equation of Figure (3-4) is: 
3 2 28 

S + (4Pp)S + (IDrry)S RP CU 
ог 


S? 4 (44a4)S^ + (4a+KB)S +Ka = 0 


where «-P and ЙВту. 


44 


BEEN B K = " a - _ 
Here G - (су "та , SO eg OF е; οτι tq. U.-0. 
ѕ +45 


КОО сасіоп оғ едцаёіопѕ (3-20) уіе1аѕ: 





2 2 


В, = —K+o = -K+1600 B, = -4% + 0) U, = 1440 
C1 = 0 Ορ = wK = -40K 
D. = и - oU, = -57600 D, = аа), - о?и, - 6400 
From equations(3-21) are obtained: 
m 57600 
-K*1600 
ЕЕ 1440:-51600)-(-К11600)(6400) (3-28) 


-40К(-К+ 1600) 


For any value of K, equations (3-28) will produce values of 
Ша апа beta that provide characteristic roots at zeta-0.5 
and omega=40. However, only certain ranges of K will meet 
the steady-state error and dominance requirements. To satisfy 
emer ror consideration it is necessary that K>1000. Since 


f (s) is of order one, equations (3-25) apply and: 


E 
p 
E ou oun, 
О 2 ВЕ ὦ 2 να 
w7 Ша 
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From the characteristic equation it is seen that 


ado = 4+0 

= 5 15 
а, 4а В 
ao = Ка 


The real part of the specified roots is š%@=20. Arbitrar ma 


choosing a dominance factor of five, the dominance criterion 


becomes: fo > 5ёш = 100. To satisfy this requirement, the 
a 

Simplest form of fg will be chosen, namely f=. It is 
б) 

i K 57600K Ш 1 36K 

then seen that 15 — ешт “ел на сч Пеи 1600=K E 


This requires that K>1176.5. Since K>1176.5 also satisfies 
the error specification, a value of K-1180 is arbitrarily 
chosen. Using this value of K, it is found that: 


а=137, 8-4.3 and ғ 2101. As a check, from the expression 


f = 
10729 


f. = (4+137)-2(0.5)(40) = ΤΟΙ 
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ШОО е 3-7 (Fourth Order Characteristic Equation) 





-e = Де > 
—— —s Каз + He ща 


Базы аты 


2-- ---------- ee ш er rs ee ee — ere | 


Figure 3-95 


Compensate the system of Figure (3-5) using tachometer plus 
 бгатіоп feedback to obtain: 
NX CHaracteristic roots at zeta - 0.5 and omega = 2. 
ГЕ 12. 
е-- 
3. Dominance of specified roots. 
The Characteristic equation of Figure (3-5) is: 


S55116.5874(73*2)8^*(82.54 8)84254K + 0 


K K 
Е S c O 
4. 54%16.559%7357%82,55%25 
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Ву inspection 6ρ” 25, e1782.5, e5773, e4716.5, EE 0571, 
U470, and Шыт . When equations (3-9) are employed one finds: 
a Ko 3 . K-24 
4 А 2 





From the characteristic equation it is seen that ад 1, аз 


а„=7З+@, a,=82.5+8, ag729*K. Since the quotient polynomial 


f (s) is a quadratic, i.e., S“+£,S+£,=0, equations (3-27) apply. 


It would be desirable that, from a dominancy viewpoint, 


51325191 = 5. However, looking at the dominancy equations 


for this case (equations (3-27)), it is seen that one of 


the several expressions for f. is fj,-a4-2549,, which is a 


fixed constant even though the remaining expressions for {1 


involve one or more variables. Thus, f,=14.5. Noting the 
a 

most simple expression for fg in equations (3-27), ΠΠ 3, 
ω 


Now, since f,=14.5°5, a dominant situation already exists. 
However, the system's performance can be further improved 


by choosing appropriate values of zeta and omega for the 





secondary roots. From the error specification it is 
Lae a 

necessary that 22712 Or K>300. Now е. (в)-57 +14.55+-0 
ω 

2 


or £,(s)=S°+14.58+6.25+0.25K. For K=300, f,(s)= 


2 
52 +14.58+81.25. Therefore, 2£505-714.5, Ше 


4279. Then zeta=0.806. These appear to be reasonable 


values for 2 and o 


=81.25 implying 


ο since the secondary roots taken alone 


would produce less overshoot and a smaller settling time 


than the primary roots. Using this value of K one obtains 


for alpha and Беш: 
48 


-16.5, 


E a. 25, SS 


and since 9-KK, and B-EK,, τ. τππε το ре ο, τοσο while 


i t 


K =0.46. 
a 
As an added bonus of this method, all the roots of the 
characteristic equation are now known and the system's 


time response could be computed if desired. 
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IV. PROGRAM DESCRIPTION 


The goal in developing any computer aided design program 
should be twofold: (1) provide the user with a single, 
easily understandable, easily usable and comprehensive 
engineering too}, and (2) dramatize its efficiency above 
that of other currently available methods. Through the examples 
of the following chapter the second goal will be demonstrated. 
It is first desirable to reveal the methodology and internal 
structure of the parameter plane curve program--as a 


consequence it is hoped that the first goal will be affirmed. 


A. THE PROGRAM 

The parameter plane curve--generating program, or 
"program" as it will be called henceforth, consists of a 
large driving routine which includes all necessary 
calculations with which to generate the curve data, and 
several supporting subroutines (i.e., curve plotting, root 
solving, data saving, etc.). This entire package is 
included as a user-selected option within another controls 
system computer aided design package. Among other options, 
the latter CAD program includes a root-locus analysis-- 
as mentioned earlier, the usefulness of either the parameter 
plane or root locus technique for design of a controls 


system is somewhat limited, but in combination their 
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КО ООС т уепесд< 15 synergistic (Chapter У will assert the dual 
roles of the root-locus and parameter plane methods). 
Facilities available within the program are many; the 

major options are: 


1. Plotting of constant zeta curves, with zeta as a 
function of omega. 


2. Plotting of constant omega curves, with omega as a 
function of гета. 


S Plotting of constant sigma (real root) curves. 

4. Plotting of constant zeta-omega curves. 

Tabular output. 

6. Rescaling of the plots. 
Input of certain data is required to enable the program. 
These inputs include: 

l. Starting value of Qe. 

2. Decades of ω to be considered. 


3. Number (and values) of constant zeta, omega, sigma, 
and zeta-omega curves. 


4. Coefficients associated with the constant, alpha, and 
beta terms of the characteristic equation. 


Each of the basic program option areas will be described in 
appropriate detail, as well as their interaction with the 
Mout data. 
1. Constant Zeta Contours 
In practice design specifications for control 
Systems are given in terms of percent overshoot, settling 
time, error constraints, etc. A value of zeta can be 


associated with the first of these specifications--that is, 


oul 


a given percent overshoot requirement сап ре related toa 
specific value of zeta. Given a specific zeta value, the 
program calculates the alpha and beta coefficients of the 
characteristic equation by holding the zeta value constant 
and varying the value of omega. The limits within which 
omega is varied are defined by the user's choice of the 
initial oi value, and the number of decades of omega to 

be considered. 

From the nature of the mapping process, it is clear 
that when the contour of the coefficient plane passes through 
a designated point (M-point), the original mapping contour 
on the S-plane passes through a point which is a root of the 
characteristic equation. The zeta value chosen for the 
contour is then the Zeta for the root. The value of omega 
associated with the M-point is the radial distance from the 
origin of the S-plane to the root. Thus, a complex root 
is determined when the M-point lies on a constant zeta curve 
Of the parameter plane. The value of this root and its 


complex conjugate is: 
= 2 
S = -£W+jwov 1-6 


If the characteristic equation is such that several complex 
roots exist, then the parameter plane curves required to 
realize these roots must all pass through the M-point. If 


г 


the complex roots have the same zeta but different omega 
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пес, then the constant zeta curve must pass through the 
Point more than once. In fact, once any point on the 
zeta contour is defined, omega, alpha, and beta are also 
defined, and all roots of the characteristic equation are 
thus fixed. 
2. Constant Omega Contours 
Within the program, for a given value of omega, zeta 
1S varied between zero and one inclusively while omega is 
held constant. 
ul fino Ons t Zeta curves, any point on a 
constant omega contour is the omega for a complex root of the 
characteristic equation. By selecting an operating point, 
zeta is also defined whereby a pair of complex conjugate roots 
is established. Again, once any point is chosen on either 
a constant zeta or a constant omega curve, all roots of the 
characteristic equation are established. 
mee Constant Sigma Contours 
When real roots are to be evaluated, it is 


Convenient to return to the characteristic equation: 


m 
БИЕ 57 - 0 (221) 


where, again, a, represents a linear combination of constant, 


k 
alpha, and beta terms. =-с (а real number) is then an 
equation of a straight line on the parameter plane. If any 


line of constant sigma value passes through the M-point, 
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then the alpha and beta coordinates of the M-point satisfy 
the characteristic equation for a real root Located ас шш 

For program considerations, one enters a positive 
value for sigma (corresponding to a real root at -σ) and 
the constant sigma contour (straight line) is developed. 
The coordinates of any point on this curve produce a real 
root at -c. That this is a useful tool, consider that 
system specifications can be achieved by placing a pair 
complex conjugate roots of the characteristic equation at 
a specific location. To ensure dominance of this root 
pair, the real part of the complex roots so placed should 
be smaller in magnitude than that of the remaining system 
roots. Roots placed at a specified Sigma value can thus 
ensure at least one real root whose magnitude is greater 
than the real part of the intended complex conjugate pair. 

4. Constant Zeta-Omega Curves 
For a fixed value of zeta and omega a pair of complex 


conjugate roots is defined in terms of the expression: 
ва 
S = -έω zw y 1-< 


The real part of these roots is, thus, defined by the 


zeta-omega product. Note that settling time is defined as 


S I . If the €w product is known, so, too, is the 


duration of the transient response. Thus, by specifying a 


o4 


constant value for the zeta-omega product, any point on the 
EMPour generated by this value will produce a given 
settling time. 

For any of the parameter plane contours, it is desirable 
Ku certain the values of the characteristic roots for every 
few values of alpha and beta. This feature has been 
incorporated within the tabular output facility as described 
below. 

EE Tabular Output 

For each of the zeta, omega, and zeta-omega parameter 
plane curves, an arbitrary though reasonable 300 points are 
calculated with which to plot the contours. For the constant 
sigma curves, only two points are needed to define the 
required straight lines (in practice, 4 points are generated 
to ensure that the sigma contours can be plotted within the 
user-defined axes limits). Because of the bulk of data 
points so generated, tabular output is offered as an option 
(as is graphical output), and all points so generated are 
ШЕГШ Са fOr the user. In addition, it is worthwhile to 
calculate system roots for given values of alpha and beta. 
However, computation of roots for each alpha and beta pair 
would cost unnecessary computer time and will likely tax the 


user's patience with the bulk of output so generated. Thus 


99 


the characteristic roots are generated for every tenth? 
pair of alpha and beta values. 
6. Plot Rescaling 

Regardless of the plot scale selected by the user, 
the program generates the full 300 data points for each 
curve requested (4 points for constant sigma lines). When 
the graphical output option is requested, the first family 
Of curves is automatically scaled to encompass each and 
every data point. The disadvantage of this technique is 
that, because most activity for the vast majority of systems 
occurs near the physical origin (i.e., alpha=beta=0), the 
curves may at first appear within only a very small sector 
of the entire plot area, and often they are indistinguishable 
from one another. The advantages of automatic scaling for 
the first set of parameter plane curves far outweigh this 
disadvantage. First, by plotting all available data points, 
the possible limits for alpha and beta are exposed--this is 
important if very large values of alpha and/or beta are 
required to meet the design specifications. Second, for 
some systems the area of activity may not occur near the 
origin, and automatic scaling spares the user the task om 


selecting a sector and possibly missing a sector of interest. 


| Although a seemingly arbitrary choice, the generation 
of characteristic roots for every tenth alpha, beta pair 
produces a very tidy output on the commonly-used IBM-3278 
computer terminal. 
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Once the user is able to view the panoramic alpha, 
beta parameter curves, it becomes obvious which sector(s) 
are of interest. The user then has the option to rescale 
the set of curves by selecting upper and lower limits for 
the alpha and beta axes. He may continue to rescale the 
family of parameter curves as often as is desirable, and at 
any time. the autoscaling option may be recalled. 

The curves generated by the above program are 
sufficient to explore most control system engineering 
problems. The use of these parameter plane contours, and 
their interaction with one another, will be evidenced in 
the next chapter. Тһе source code listing of the program 


is included as Appendix B. 


KS INSTRUCTION TO THE USER 

The parameter plane program is highly interactive--the 
user is prompted for each required input. A brief 
description of all but the most trivial input items follows. 


- Starting value of ® : For most control systems the 
initial value of ? is chosen to be zero. Because 9 
is used in the denominator of certain of the parameter 
plane equations, ὦ must be greater than zero. However, 
the user may choose © arbitrarily close to zero if 
desired. $ 


- Number of decades: For the majority of control system 
problems a suitable number of decades to be considered 
might be two or three. For higher order systems, it 
would be advisable to start with a slightly larger 
number of decades, especially if the initial ? value 
is small. For subsequent families of curves, the 
number of decades can easily be changed. 


OT 


Тһе 


4. 


Constant, alpha, and beta coefficient values: 
Characteristic coefficients are requested from the 
highest to lowest order term. By way of an example, 
a third order characteristic equation might be: 


53 + (За + а + 10)52 + а5 + (8+5) = 0 


Here, the constant coefficients would be entered in the 
following sequence: 1,10,0,5 while alpha and beta 
coefficients would be entered as 0,3,1,0 and 0,1,0,1 
respectively. 


Zeta values: By convention, values are restricted to 
between zero and one, inclusively. 


Sigma values: Positive values of sigma correspond to 
negative real roots. Since few, if any, practical 
engineering applications exist for designing a positive 
real root into a system, negative values for sigma are 
disallowed. 
Omega values: Values for constant omega curves are 
restricted in the lower limit tQ the starting s value, 
and in the upper bound by o, x10 SE ай 
Zeta-omega values: As with the constant sigma curves, 
values for constant zeta-omega contours must be greater 
than or equal to zero. 
user then has the following options: 

Review entries. 

Change any entry. 

Tabular output. 

Graphical output. 


Remember that tabular output includes 300 data points 


for all but the sigma contours. Characteristic roots are 


displayed for every tenth alpha, beta pair. Because of the 


bulk of output for this optiomi use it only when necesse jg 


If a printed copy of the tabular output is desired, type 
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"record on" before invoking the program. Upon exiting the 
pnssram. type "record off", after which the user may save 
the preceding terminal session in a listing file designated 
ша пале от Піс choosing. Simply print the listing file, 
which will include all output which has transpired on the 
terminal between the two calls to "record". 

When graphical output is requested, all curves are 
Superimposed on the same plot. The first set of curves is 
produced with an autoscaling feature, which plots all 
points calculated (the range of points depends on your 
choice of initial Ee number of decades, zeta values, omega 
values, etc.). For most characteristic equations only the 
first quadrant (i.e., positive alpha and beta values) will 
be of interest, since negative values usually (but not 
necessarily) imply negative characteristic coefficients and, 
thus, positive roots leading to system instability. Тһе 
nature of the first (autoscaled) set of curves will reveal 
the actual areas of interest for subsequent plots for the 
same system. 

Finally, after each family of curves is plotted, the 
user has six additional options: 

1. New problem. 
2. same problem. 
Ec Root finder. 


4. Save problem. 
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5. Create "DISSPPA meta irers 

6. Return to main menu. 
Items 1, 4, and 6 are self-explanatory. The remaining 
options deserve some additional explanation. 


- Option 2: This option can be used to re-enter the 
problem at a point prior to actual plotting. Then, 
specific input values can be added or modified, 
tabular and/or graphical output can be requested, and 
entries can be reviewed. It is within this option that 
the graph coordinate axes can be rescaled. If user- 
defined scaling is desired, the minimum and maximum. 
values for the axes are requested. 


- Option 3: Although within the tabular output feature a 
set of characteristic roots is produced for every tenth 
pair of alpha and beta values, the bulk of output using 
that option may prove excessive for some applications. 
Here, the user has the option of choosing specific 
values of alpha and beta (e.g. extracted from the 
family of parameter plane curvés) and obtaining the 
system roots. 


- Option 5: The program provides a choice from among 
four graphic output devices. Usually the user will 
nominate the TEK618 graphicsterminal due to its 
relatively high quality plot resolution. Once the 
user has produced a parameter plane plot to his liking, 
he may wish a final plot of very high resolution. By 
selecting this option, the plot is stored as a DISSPIA 
metafile, and the program is terminated (termination 
of the program is necessary at this point due to an 
anomaly of the DISSPLA graphics package). Simply type 
"DISSPOP" and follow the simple instructions, choosing 
the default options as they are presented. Within the 
"DISSPOP" routine, any of several output devices can be 
called, including the high resolution Versatee plotren 
and the 3800 laser printer. 
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V. PARAMETER PLANE CURVES-GRAPHICAL METHOD 


A. GRAPHICAL SOLUTION 

The algebraic solutions discussed in Chapter III have the 
disadvantage that a fixed value of zeta and omega must first 
NE chosen to compute the alpha and beta terms. In some 
instances it is possible to modify the remainder polynomial 
SO as to ensure that the specified roots are dominant. 
However, it is not always possible to guarantee that roots 
placed at a specified location can be made dominant. Thus, 
an exhaustive trial-and-error procedure may be required to 
achieve the best values for the various parameters.  Trial- 
and-error may also be required in the design of cascade 
compensators where a specific root location may require 
parameter values that are not physically realizable. In 
these cases, the calculation must be repeated in terms of 
EE. modrtred Specifications; possibly a different means 
Of compensation must be used. 

Bu Ота ои trial-and-error analytical approach, one 
can employ the graphical solution. Once a family of curves 
is generated by the program one can, by choosing an M-point 
in the parameter plane, obtain from the curves the n roots 
of the nth order characteristic equation. The trial-and-error 
procedure can then be done visually to reveal an operating 


point which best meets the given specifications. 
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Example 5.1 (Àn Attitude Control System for Large Launch 
Vehicles) 


RIGID BODY 
DYNAMICS 





CONTROL CONTROL 
COMPUTATION LAG 


STABLE PLATFORM 





ACCELEROMETERS 





` Figure 5-1 


| 
| 

GUIDANCE] ` OTHER 
LAW SENSORS 


Figure 5-1 shows the mechanization for a control system for a 


large launch vehicle; Figure 5-2 shows the equivalent block 


| 

diagram 
| | 

| | 

| 





where the equivalent G(s) is: 


(agS? « K,S « K,X(S^-C) 
- T 2 
ο τω τν o 
ο δι ο) 
(Ú (Ú 1 


5 ero 4 | 3 2 
= — 2a ха 
"z*àü ο. (1*ag)S + (a, C*K, )8 * (K, а „С)5 
е 
е 
- СК. 5 - CK, = 0 
мһеге Δρ. 81; Κι Κο - control system gains 
A = damping ratio. of control servo 
E = - natural frequency of control servo 


Gains K; and Ko are chosen as the system parameters and 6, 
respectively, to be portrayed in the parameter plane. One 
must then find suitable values for these parameters to yield 
a desired stability margin and a satisfactory transient 
response. For a typical choice of system parameters (such 


as those describing Saturn V), it is assumed that: 


о О 7 

е 
P. = 4.71 Hz = 29.594 radians 
а, 75 
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Then the всувТеп в спагасгектеріс еппагитивесошесі 


0.001156 + 0.04858? E ο.55΄ ш (о.т+а)5° πο ο > T 8)s* 


_ 0,76 - 0595. 0 


Various values of ag and а, were used to deduce their effect 
on the stability regions, which is indicated in Figure 3-3. 
The numbers of stable and unstable roots, respectively, are 
portrayed in parentheses for each region of the parameter 
plane. 

The analysis procedes as follows: (һе €=0 contour осте и 


sents boundaries of stability (or relative stability when 


€>0) associated with pairs of complex conjugate roots. The 


S-0 (sigma-0) curve represents real root stability boundaries. 


The region of stability is thus that area bounded by these 
two contours. See Figures 5-4a and 5-4b for a magnified 
view of this area. Any negative alpha-beta pair from within 
this region will exhibit six stable (i.e., all within left- 
half S-plane) roots and system stability will be assured. 
Note that from the form of the characteristic equation, 

both alpha and beta must be negative to obtain a stable 
System. To illustrate, let us select an arbitrary operating 


point constrained to lie within the lower loop of Figure 5-4b 
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Parameter Plane Curves for 
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Figure 5-4a 
G Parameter Plane Curves for 2 
0.0011S° + 0.0485S° + 0.5S* + (0.7+a)S% + (0.35+8)S 
- 0.785 - 0.70 = 0 


——— = а 
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ATTITUDE CONTROL SYSTEM 
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Figure 5-4b 
Parameter Plane Curves for 3 2 
0.001186 + 0.0485S° + 0.584 + (0.7+a)S~ + (0.35+8)S 
20m as Ong c= = 
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and also satisfying the requirement that both alpha and beta 
be negative. Superimposed on Figure 5-4b is the constant 
€=0.5 contour, upon which our M-point might be chosen. If 
we select, say, %=-0.15 and §=-0.02 (corresponding to Е +0.5 
and 920.5), the system roots are: 

-0.26 + 30.45 

-dominant roots 

-0.26 - 10.45 

-0. Заа 50 2 

-0:32-- 10.12 

-14.6 + 10 


-20 EO 


Coincidentally, the roots associated with our choice of zeta 
and omega are seen to be dominant. The actual choice of an 
Operating point may depend on other criteria not discussed 


here. 


Since Ky and K, (¿ and 8 , respectively) are functions 
ONE 5 ag» 84» С, ш, and £,, they тау be determined for 
various instances of flight by plotting several constant 
zeta and constant omega curves. Actually, K; and Ко vary SO 
little within the range of values used for C, for specified 
values of 5 and 9, that it becomes possible to choose constant 
values for K; and Ко. 

This has been a relatively simplistic treatment of a 


complicated control system problem, but it demonstrates the 


power of the parameter plane graphical technique. Knowing 
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nothing more than the system's characteristic equation, the 

static and dynamic stability boundaries may be obtained to 

define the area of overall stability. Of course, other 

System constraints may exist to further limit this area. 

As a note, the root-finding option available within the 

program was used to confirm the numbers of stable and 

unstable roots lying within each region of Figures 5-3 and 5-4. 
Example 5-2 (Alternator Voltage Regulator) 

In designing a voltage regulator of the type shown in Figure 


9-9, we must find values for К K and K, that provide 


Q? ДА 


stability and good transient performance. 


е Сем ALTERNATOR 






2 
KS *K,S*Kg 
(TS 1)(T,S 1) 








Figure 5-5 


The characteristic equation of the system, including alternator, 


tie-line, etc. is (see reference 8 for details): 


1850955?» 0.1325S? + 1.728” + К с 5: + (K +9.1)S 


+ Ку -2.9 = 0 
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We will consider Жы as fixed and o and 8 as variable para- 
meters corresponding to K,*7.55 and K,*9.1, respectively. 
Curves of а versus P are plotted in Figure 5-6 with a, as 
the variable parameter (for this systen, κι апа K, are 

known to be functions of ως alone). Since the program plots 
constant zeta curves as a function of varying omega, these 
curves were selected as logical candidates to study the 
problem. 

The beta axis corresponds to the zero value of damping 
constant (sigma) since below the zeta-1.0 curve, the real 
roots (Thaler and Brown 1960) are the negative slopes of the 
tangents drawn from the point in question to the £-1] curve. 
The machine is then stable for any (a, 8) pair between the 
650 and S-0 curves. Тһе greatest stability of the machine 
is then possible when both zeta and sigma are largest. 
Further, the best stability can be expected in the region 


bounded by the £-0.3 and £-1.0 contours (Kabriel 1967). 


Similar stability limits о? а апа 8 can be investigated by 


taking Ko as 10, 20, 30 and so on. 
Consider now K,=0 as fixed. Then a and 8 will represent 
the variable parameters K,+7.55 and K,-2.5, respectively. 


2 О 


The characteristic equation then becomes: 


0.00958? «4 0.1325S* « 1.728? + aS? + 9.1S + 8 = O 
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Again, Ko and Ko are known to be functions of ae alone. Here, 


however, € can be solved explicitly [Ref. 8] and one obtains 


three straight-line equations: 


B = O 
5.421% -β = 3.894 
175.63a -6 - 4085.0 


or in terms of system parameters: 


Ko 


5.421(K,+7.55)-K 


= 


3.894 - 2.5 


175.63(K,*7.55)-Kg 4085.0 - 2.95 
These are plotted in Figure 5-7. The triangular region 
bounded by these three lines represents a stable region. The 
values о (а, В) within the triangle and hence corresponding 
values of Ko and Ko can be predicted for stable operation. 
The procedure can be repeated for further investigation by 
taking Ку=15, 30, 45, etc. 

The analytical parameter plane technique can be used to 


determine the stability limits of Ko and K, by choosing several 


1 


constant values of Ко: But since K, has to be selected 
arbitrarily for this purpose, this method becomes cumbersome 
and time-consuming. The method presented here, on the other 


hand, can be used to determine the stable range of either 
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Figure 5-7 
5 Parameter Plane Curves for 
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9 


Ko and Ко ог κι апа K, fixing the third parameter. Over επ. 
above the ability to predict stable operation, the method 
povides a direct measure of the damping at and around a 
chosen operating point. 


Example 5-3 (Lead Compensator) 


Consider the system: 






1 
$(5+1)(5+4) 


The characteristic equation is: 


st + (5+р)83 + (4*5P)S* + (4P+K)S + KZ = 0 


We choose to cancel the pole at S=-1 with the zero; thus 


Z=1.0 and the characteristic equation becomes: 


oh Sip) Snes" Sp Soren = ae 


Let P =a and K = В. Then the parameter plane curves are as 


shown in Figure 5.8, Form the coefficients of өне сһагас сой 


istic equation to remain positive (and thus ensure stability), 


it is convenient to consider only positive values for alpha 


” 


and beta. By inspection we choose §=0.9 and 0152.09 аз а 
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Figure 5-8 
Parameter Plane Curves for 
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Figure 5-9 
Parameter Plane Curves for 


+ νο. + (44+50)8° + (4048 )g +8 
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О 


good operating point, for which %® = 4.0 апа В -24.0. The 

corresponding vicinity of Figure 5-8 is re-scaled in Figure 

0-9. Again, using the root-finding option, the roots 

r Sc l ¿od withubhis,. (S. gs) pair are shown to be dominant. 
Example 5-4 (Lag Compensator) 

If we are especially concerned with steady-state accuracy 

for a ramp input, it may be advisable to design a lag 

compensator. The parameter plane permits us to consider 

steady-state accuracy while designing the transient response. 

If our system is the same as that considered for the lead 


compensator, the characteristic equation remains: 


3m + (5+P)s° + (4155γ5΄ + (4P+K)S +KZ = 0 


But now the error coefficient is K = 54 . Having three 
unknown parameters, K, Z, and P, two must be selected (or 

some combination of two) to be Q and P while a numerical value 
is assigned to the third. Once this choice has been made 

the parameter plane curves can be calculated and plotted, and 


the loci of constant K can be superimposed. Let 7-0.1, 


Р-0 апа Ке8, The characteristic equation becomes: 


54 T (5+ 5 ri (4*50)g* + (40+-8)S + 0.18 = 0 


Parameter plane curves are shown in Figures 5-10 and 5-11. 


В 
— - 0.1, 0.2, 0.3, etc. may be superimposed. 


If we select K = 0.15 and (similarly to the lead compensator 





Lines of K - 
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Figure 5-10 
A arameter Plane Curves for 
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GO 


example) €=0.5 and в =2.0, then 0=2.14 үһі1е Е-12.89, This 
produces dominant roots atr Pm 4 11,7. 

Let us reconsider the systems for which a lead compensator 
and a lag compensator were designed. The characteristic 
equation as well as the error coefficient each contain three 
unknowns (if we assume a numerical value for К). Че сап 
imbed the error coefficient in the characteristic equation 


by direct substitution, thereby eliminating one of the 





unknowns. Let us eliminate the gain parameter K - note that 
4PK 
К = 7 r 


Returning to the characteristic equation: 


4 4PK 


S + (5+Р)5° + (ДЕ + (4р+ = SS + АРК = ¢ 





Letting P = а апа 5 - В, the characteristic equation becomes: 


S* « (541)8? « (4«50)S^ « (4u*4K B)S * 4K a = O 
If à constant value is now chosen for K Say К = 2.0, the 


characteristic equation becomes: 


st 4 (54a )S° + (4450 )8° « (dn488)S 4 80 = 0 


for which the parameter plane contours are first shown in 
Figure 5-12 and are further magnified in Figure 5-13. Now 
when any operating point is chosen on the parameter plane 
curve(s), the selected (9,Pb) pair generates K 72.0. Of 


course, this procedure can be repeated for any choice of Ка: 
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A Parameter Plane Curves for 
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Figure 5-13 
A Parameter Plane Curves for 
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On the parameter plane curves let us choose 5=0.5 and 


22.9 LOr Which: 





m 6,00 = Pp 

Е ae 
ЕЕ 1.50 - 2 
α 
= 1359 = Z 
К + 83 = 36 

: τν 
As a check, κ. s up 2.0. The roots associated with the 


selected zeta and omega values are shown to be dominant 
uSing the root-finding algorithm and are: 
moO oF 1.73 
- Dominant Pair 
ee OO = 11.72 
63 Ғ 10 


feat + 340 


Example 5-5 


AS a final engineering example, consider the system below: 


+ K 
| е — ~|S(S+10)(S?+5S+100) 
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For this system, the following specifications are to be meses 
1. Set K to the stability апып. 


2. Place a dominant root pair within the following region: 
0.4<5<0.7, ала 2<ш<6. 


3. Both tachometer and acceleration feedback may be used, 
but if possible choose only one. 


From the figure the uncompensated system's open loop transfer 


functions: 


K 


GH = sl EURO s 
S(S+10)(S +5S+100) 


From which, when expanded, the characteristic equation 


becomes: 


54 s 158° + 1508“ + 1000S + K = 0 


To determine the value of K at the stability limit the Routh 


array is employed: 


1 150 K 

15 1000 0 

1250 15K 0 
1.25х109-225К 0 0 
15К 0 0 
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Шого (пе stability limit is seems to be K=5555.5. If both 
tachometer and acceleration feedback are used the compensated 


system's characteristic equation becomes: 


Si E 155° τ ТООЛО S ο 0005555725 ο NDDOD.O-0 


where aSK > B=K > and K has been set to the stability limit. 
Parameter plane curves for this system are shown in Figure 
20-14. 

From these curves, the following analysis can be made. 
The origin of the parameter plane corresponds to the roots 
of the uncompensated system. Since the 5-0 curve passes 
iieougi the origin, two roots*are located on the jw axis of 
the S-plane as was to be expected from the Routh array. The 
remaining two roots are also complex and correspond to &=0.8 
and w=5.0. It is important to note that when an operating 
point involves two different pairs of complex roots, then the 
curves for two different values of omega and two different 
values of zeta must pass through the point. To determine 
which value of omega corresponds to which value of Zeta, it 
becomes necessary to refer to the program's tabular data 
ОШКОШ, which is not included here due to lack of space. 

With K =O, the effect of tachometer feedback alone 
corresponds to movement of the M point along the 6 axis. 
In Figure 5-14 the unstable region is determined by an 


K Section of the way the constant zeta curves tend 
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ο л псгеалсесз. Since the B axis is always in the 
unstable region, it is concluded that tachometer feedback 
alone cannot stabilize the system. 

The effect of acceleration feedback alone can be observed 
EBXravelding along the a-axis of Figure 5-14. If K is 
varied between 0.01 and 0.06, the system will exhibit two 
pairs of complex roots with the following ranges of values 


for zeta: 
шие 0 5 апа 0.25<5<0.32 


If both tachometer and acceleration feedback are used, 
it is seen that tachometer feedback will in general cause 
the zeta of one pair of roots to increase while the zeta 
of the other pair decreases. Acceleration feedback alone 
would appear to be the better choice. 

From the set of curves it is determined that with K, -0 
and K 70.012, four complex roots are located with associated 
ΜΙ ас 5=0.45, w=4.0, and £20.32, w=13.0. Since 
(0.45)(4)218:*(0.32)(13)24.15, it is apparent that the roots 
at £20.45 and 924.0 are dominant. The specifications have 


been met and the problem is solved. 


B. MISCELLANEOUS ASPECTS OF THE PARAMETER PLANE 
It can be demonstrated that constant zeta parameter 
plane curves of order two through five originate at a point 


where alpha - = and beta = = where M, N, and K are 
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determined by the zero and first power coefficients only. 
Intuition can be used to conclude that constant zeta curves 
of any order originate at this. Common point which S 
determined only by the zero and first power coefficients. 

An exception is when K=0. In this case the origin point 
depends on higher order coefficients and its location ІІІ 
be obvious given a specific problem. If K is not zero, 

the origin is independent of the order of the characteristic 
equation. 

Inspection of the expressions for alpha and beta indicates 
that the shape of the constant zeta curves as omega becomes 
larger is primarily determined by the coefficients of higher 
power, and in general the curves become more complex and 
less well behaved as the order of the characteristic 
equation increases. For a given characteristic equation, 
an increase in complexity can be observed as alpha and beta 
appear infmore coerficients.: 

А11 constant zeta curves tend to plus or minus infinity. 
The relative magnitudes of the coefficients determine whether 
the limit is plus or minus infinity. It is therefore 
necessary to choose a frequency range of interest before 
plotting the curves, спивя limibime the analysis to one 
"window" of the infinite plane. 

Since no stability criteria, either relative or 


Ж 


absolute, has been established for the parameter plane, it 
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is necessary to base the stability analysis on observing 
which way the curves tend as omega and zeta are varied. 

For this reason it is worthwhile to plot curves for as many 
values of zeta, omega, Sigma, and if desired, zeta-omega, 


as are necessary to ascertain the pattern. 


Go 


ҮІ. CONCLUSIONS 


Parameter plane techniques have been applied to the 
compensation of linear control systems. General equations 
have been derived for the cases of feedback, cascade, and 
combination feedback-cascade compensation, to enable one to 
place a pair of complex conjugate roots at a specific 
location in the S-plane, while simultaneously satisfying 
the steady-state accuracy requirements. A dominancy 
technique has been introduced whereby once a pair of complex 
roots is fixed, the remaining roots of the characteristic 
equation can be manipulated to ensure that the specified 
roots are dominant. y 

The impetus for development of a parameter plane program 
was to provide the user with a quick, simple means of 
obtaining the information available in the analytical 
solution to control system compensation, while avoiding the 
painstaking labor of trial-and-error analysis inherent in 
that technique. Several, practical engineering examples 
have been presented to demonstrate the superiority of the 
graphical technique. To date, no other package is known to 
offer the fully interactive and comprehensive capabilities 


of the parameter plane program. 
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By itself the program allows one to design a control 
system compensation model for most systems. However, for 
some lightly damped systems containing mechanical resonances, 
the amount by which zeta or omega are incremented in the 
parameter plane equations may be so large as to not detect 
the resonance peaks. This information would be available from 
either a root-locus or Bode analysis. For still other 
systems, one might be interested in the way the roots of the 
characteristic equation extend from the open loop poles and/or 
zeros. Since the parameter plane equations are calculated 
uSing only the characteristic equation, no Knowledge of open 
loop poles or zeros is available. Again, a root-locus method 
would reveal this information. Incorporated into one 
comprehensive package which includes Bode and root-locus 
analyses, the program provides the capability to investigate 
the entire gamut of linear control system architecture. 

A basis for further investigation involves plotting the 
parameter plane contours for systems that are non-linear 
in the alpha and beta terms--i.e., those systems which contain 
alpha-beta product terms. Although the recursion technique 
used in this text has distinct advantages over the matrix 
approach for the linear case, as pointed out earlier, the 
matrix technique would be the method of choice for the non- 


linear case. 
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APPENDIX B 
PARAMETER PLANE PROGRAM 


SUBROUTINE LPARAM 

DIMENSION A(350), B(350), 
AG(350), BG(350), AJ(100), BJ(100), CJ(100), 
ZETAC100), SIGMAC100), WC100), ZWC100) 


INITIAL ASSIGNMENTS 

CHARACTERX4 5СНАК/"5 '/,YES/'Y '/,NOO/'N '/,BLANK/' "И, 
ENCHAR/'E '/,MWNCHAR/ 'WN'/,NDCHAR/ 'ND'/,NOCHAR/''NO'7, 
AJCHAR/'AJ'/,BJCHAR/'BJ'/,CJCHAR/'CJ'/,NSCHAR/ 'NS''7, 
NZCHAR/ 'NZ'/,ZWCHAR/' ZW'/ , NHCHAR/ 'NW'/, PRCHAR/ 'PR'/, 
NCCHAR/ 'NC'/,TICHAR/'TI'/ 

CHARACTERX4 TABLE, GRAPH, CHANGE, REPLY, OPT, LABEL(C9) 

COMMON /SAVE/ LABEL, WN, ND, N02, NC, CJ, AJ, BJ, 


NZ, ZETA, NS, SIGMA, NW, W, NZW, СИ, 
XMIN, XMAX, YMIN, YMAX 
DATA ENTRY FROM FILE OR CONSOLE? 
MINMAX = 1 
GRD = 0. 


CHANGE = BLANK 
CALL EXCMSC'CLRSCRN') 
WRITEC6,500) 
CALL READC CREPLY) 
ШЕ CREPLY .NE. 'D') GOTO 101 
CALL GETIT 
MINMAX = 0 
GO TO 200 
CONTINUE 


GET CURVE TITLE 
CONTINUE 
CALL EXCMSC'CLRSCRN') 
WRITEC6,501) 
CALL READL (LABEL) 
CALL ASTER (LABEL,LABEL) 
IF ( CHANGE .EQ. TICHAR ) GO TO 200 


GET STARTING VALUE OF WN 

CONTINUE 
МКІТЕ(6,502) 
CALL READR (ИМ) 

IF (ΝΝ) 104,104,105 

WRITEC6,503) 

GO TO 103 
IF € CHANGE .EQ. WNCHAR) GOTO 200 


GET THE NUMBER OF DECADES CONSIDERED 
CONTINUE 
WRITEC6,504) 
CALL READI (ND) 
IF ( CHANGE .EQ. NDCHAR ) GOTO 200 


GET THE ORDER OF THE CHARACTERISTIC EQN 
CONTINUE 
WRITEC6,505) 
CALL READI (N02) 
NC = М02%1 
IF (CHANGE .Е0. NOCHAR ) GOTO 200 
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O O 


ОО 


OM 


ОООО 


108 


109 


110 


111 


115 


116 


117 


118 


119 


GET THE NUMBER OF CONSTANT ZETA CURVES 
CONTINUE 
CALL EXCMS (C'CLRSCRN') 
МКІТЕ(6,506) 
CALL READI (NZ) 
IF (NZ .LT. 1) GOTO 110 


GET THE VALUES OF ZETA 

WRITEC6,507) 
DO 110 I = 1,NZ 

WRITEC6,508)I 

CALL READR CZETACIO) 

ТЕ (ФЕТА(Т) .LT. 0. „ОК. ФЕТА(Т) „СТ. 1.) WRITEC6,509) 

ТЕ (ФЕТА(Т) .LT. O. „ОК. ZETACIO σσ οὐ τοι 
CONTINUE 
ТЕ ( CHANGE. EQ. NZCHAR ) GOTO 200 


GET THE NUMBER OF CONSTANT SIGMA CURVES 
CONTINUE 
CALL EXCMS C'CLRSCRN') 
WRITEC6,510) e 
CALL READI (NS) 
IF (NS .LT. 1) GOTO 113 


GET THE VALUES OF SIGMA 


DO 113 I = 1,NS 
WRITEC6,511) I 
CALL READR CSIGMACI)) 
IF CSIGMACI) «ΙΤ. 0.) WRITE (6,512) 
IF CSIGMACID .LT. 0.5 GO TO 11? 
CONTINUE 


IF € CHANGE .EQ. NSCHAR ) GOTO 200 


GET THE NUMBER OF CONSTANT ИМ CURVES 
CONTINUE 
CALL EXCMS C'CLRSCRN'!) 
WRITEC6,515) 
CALL READI (МИ) 
IF (NW .LT. 1) GOTO 116 


GET THE WN VALUES 
WNMAX = WNX10xXxND 
DO 116 I = 1,NW 
WRITEC6,514) I 
CALL READR (И(Т) ) 
IF (NCI) .LT. ΗΝ .OR. WCI) .GT. WNMAX) WRITE (6,515) ММ, ИММАХ 
IF (NCI) .LT. WN .OR. WCI) .GT. WNMAX) GO TO 115 
CONTINUE 
IF € CHANGE .EQ. NWCHAR ) GOTO 200 


GET THE NUMBER OF CONSTANT ZETAXWN CURVES 
CONTINUE 
CALL EXCMS C'CLRSCRN!) 
WRITEC6,516) 
CALL READI (NZW) 
IF (NZW .LT. 1) GOTO 119 


GET THE ZXMWN VALUES 
DO 119 I < 1, МЕН 
WRITEC6,517) I 
CALL READR (ZWCI)) 
IF CZNCID) .LE. 0.) WRITE (6,518) 
IF CZNCID .LE. 0.) GO TO 118 
CONTINUE 


“ДЕ © CHANGE „Еф. ZWCHAR ) GOTO 200 
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O O ОО 


OM 


ОООООООООООООООООООООООООООООООООООООО 


120 


121 


122 


125 


124 


125 


GET CONSTANT COEFFICIENTS 
CONTINUE 
CALL EXCMS С "СІ КСКМ") 
МКІТЕ(6,519) 
DO 121 I = NC,1,-1 
II = I-l 
МКІТЕ(6,520) ІІ,І 
CALL READR (CJCI)) 
CONTINUE 
IF € CHANGE .EQ. CJCHAR ) GOTO 200 


GET ALPHA COEFFICIENTS 
CONTINUE 
CALL EXCMS C'CLRSCRN') 

МКІТЕ(6,521) 
DO 125 I = NC,1,-1 

= ТЕТ 

ИКТТЕС6, 522) 11, 1 

CALL READR САЈСІ)) 
CONTINUE 
IF (CHANGE .EQ. AJCHAR ) GOTO 200 


GET BETA COEFFICIENTS 
CONTINUE 
CALL EXCMS C'CLRSCRN') 
WRITEC6, 523) 
DO 125 I = NC,1,-1 
II = рет 
МКІТЕ(6,524) ІІ, 
CALL READR CBJCI 
CONTINUE 


I 
)) 


25 


С 
С 


G> O, O, O, O, G; G; G; GO; O, O O, O O O O O O O 


200 


201 
202 


203 
204 


205 
206 


207 
208 


209 


REVIEW ENTRIES 

CONTINUE 
CALL EXCMS C'CLRSCRN') 

WRITEC6,525) 

CALL READC CREPLY) 

IF CREPLY .NE. YES) GOTO 209 
CALL EXCMS C'CLRSCRN'!) 
WRITEC6,526) 

WRITEC6,527) CLABELCI), I=l1, 9) 

ККІТЕС6, 528) 

ККІТЕС6, 529) МЮ, №02, №, М5, МИ, МАМ 

WRITEC6,530) WN 

WRITEC6,531) 

IF (NZ .LT. 1) GOTO 201 
WRITEC6,532) CZETACM) ,M=1,NZ) 
GOTO 202 
WNRITECÓ,533) 

CONTINUE 

WRITEC6,534) 

IF (NS .LT. 1) GOTO 203 
WRITEC6,532) С5ІСМАСМ) , М=1,М№5) 
GOTO 204 
WRITEC6, 533) 

CONTINUE 

WRITEC6,535) 

IF CNW .LT. 1) GOTO 205 
WRITEC6,532) CWCM),M=1,NW) 
GOTO 206 
WRITEC6,533) 

CONTINUE 

WRITEC6, 536) 

IF (NZW .LT. 1) GOTO 207 
WRITEC6,532) (ZWCM),M=1,NZW) 
GOTO 208 
WRITEC6,533) 

CONTINUE 

WRITE (6,576) 

CALL EXCMSC'CLRSCRN'!) 

WRITEC6, 537) 

WRITEC6,532) CCJCN),N=NC,1,-1) 

МЕІТЕ(6,538) 

WRITEC6,5352) CAJCND,NzNC,1,-1) 

WRITEC6,539) 

WRITEC6,532) CBJCND,NzNC,1,-1) 

WRITEC6, 540) 

WRITEC6,541) XMIN, XMAX, YMIN, YMAX 

CONTINUE 

WRITEC6, 542) 

CALL READC (CHANGE) 

IF (CHANGE .NE. YES) GOTO 210 
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O; O, O, O) O) G; O; O) G; O) O; O; O; O; G; O) O; G; O; O; O; O; O; O) O O) O; O) O; O) O G; O) O, O O, GO O 


210 


CALL EXCMS (C'CLRSCRN') 


WRITEC6, 543) 
CALL READC ( 


IF 
IF 


WRITEC6, 544) 


(CHANGE 


(CHANGE . 
(CHANGE . 


( CHANGE 
( CHANGE 


(CHANGE . 


(CHANGE . 
(CHANGE . 
(CHANGE . 


( CHANGE 


(CHANGE . 


(CHANGE 


(CHANGE . 
(CHANGE . 
GOTO 209 
CONTINUE 


CHAN 
EQ. 


GE) 

TICHAR) 
WNCHAR) 
NDCHAR) 
NOCHAR) 
NZCHAR) 
NSCHAR) 
NHCHAR) 
ZWCHAR) 
CJCHAR) 
AJCHAR) 
BJCHAR) 
ENCHAR) 
PRCHAR) 
NCCHAR) 


CALL READC (TABLE) 


МКІТЕ(6,545) 


CALL READC (GRAPH) 


ІҒ (М2) 215,215,212 
- (10.XXND)*Xx(1.7299.) 


G 
CONTINUE 
2 = 1.0Е-60 


CHANGE ROUTINE 


GOTO 
GOTO 
GOTO 
GOTO 
GOTO 
GOTO 
GOTO 
GOTO 
GOTO 
GOTO 
GOTO 
GOTO 
GOTO 
GOTO 
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ооо O O O 


CONSTANT ZETA PLOTS 
IF (NZ) 309,309,300 
300 CALL EXCMS C'CLRSCRN') 
IF (TABLE .EQ. YES .OR. GRAPH .EQ. YES) WRITE (6,546) 


DO 308 M=1,NZ 


IF (TABLE .EQ. YES) WRITE (6,547) 
J = 0 
К = 0. 
ИМА = WN 
DO 306 L=1,300 
Dl - 0.0 
D2 - 0.0 
С1 = 0.0 
С2 = 0.0 
В1 = 0.0 
В2 = 0.0 
DO 303 N=1,NC 
K = N-1 
IF (K) 302,301,302 
501 U = 0.0 
Ul = -1.0 
302 U2 = 2.0XZETACM) XU-U1 
D1 = (-1.00XXKXCJCNOXHNAXXKXUT-*DI 
D2 = (-1.0)XXKXCJ(N)XWNAXXKXU+D2 
СІ - (-1.00XXKXBJCNOXWHNAXXKXUI*CI 
C2 = (-1.0)0XXKXBJCNO XHNAXXKXU-*C?2 
Bl = (-1.00XXKXAJCNOXHNAXXKXULI* B1 
B2 » (-1.00XXKXAJCNOXHNAXXKXU- B2 
Ul = U 
U = U 
303 CONTINUE 
ТЕ CABSCB1XC2-B2XC120-2) 306,306,304 
304 J = Jtl 
R = R*1. 
ACJ) = (C1xD2-C2xD1)Z(B1xC2-B2xC1) 
BCJ) = CB2*XD1-B1*¥D2)/CB1XC2-B2xC1 ) 


IF (TABLE .NE. YES) GO TO 306 
WRITEC6,548) ACJ), BCJ), ИМА, ZETACM) 
IF (R410. - 410) 306,305,306 
305 CALL ROOTS (ACJ), BCJ), AJ, BJ, CJ, N02) 
CALL EXCMSC'CLRSCRN') 
WRITE (6,547) 
506 WNA = GXWNA 


CALL EXCMSC'CLRSCRN') 
IF (J .GT. 0) GOTO 307 
МКІТЕ(6,549) 
GOTO 308 
307 IF (GRAPH .EQ. YES) CALL PLOTDCA,B,J,.FALSE., 
LABEL, 'ALPHAS','BETAS', 
МІММАХ," 255",2ЕТА(М), 
XMIN,XMAX,YMIN,YMAX,GRD) 
IF (GRAPH .EQ. YES) GRD - GRD*1. 
308 CONTINUE 


309 CONTINUE 


>< >< >< 
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С 
С 
С 


С 


Оооо O O O 


325 


326 


327 


328 


329 


330 


351 


352 
555 


334 
335 
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IF (М5) 
CALL EXC 
IF ¢ 


XINC 
YINC 


DO 


CO 
CONTIN 


CONSTANT SIGMA PLOTS 
33573397329 
MSC'CLRSCRN') 
TABLE .EQ. YES .OR. GRAPH .EQ. YES) WRITE (6,550) 


- (ХМАХ - ХМІМ)/299. 
- (ҮМАХ - ҮМІМ)/299. 


334 M=1,NS5 
= XMIN 
YPT = YMIN 
CJC1) 
BJC1) 
AJC1) 


O 
O 
пни 


DO 326 N-2,NC 
K - N-1 
DUMMY5zSIGMA(CMOXxK 


DD = (-1.00 XXKXCJCNOXDUMMY5-*DD 

CC = (-1.00 XXKXBJCNOXDUMMY5-*CC 

BB - (-1.0)0XXKXAJCNOXDUMMY5-*BB 
CONTINUE 


IF (CC .EQ. 0. .AND. ВВ .Е0. 0.) 60Т0 334 
IF (CC) 327,327,330 
DO 329 L=1,300 
= J+1 
= R+1. 
(J) = XPT 
(J) = (-BB%xA(J)-DD)ZCC 

IF (TABLE .NE. YES) GOTO 329 

МКІТЕ(6,552) А(,/), В(,/), 5ІСМА(М) 

IF (R710. - J710) 329,328,329 ` 

CALL ROOTS (ACJ), BCJ), AJ, BJ, CJ, N02) 

CALL. EXCMS C'CLRSCRN') 

WRITEC6,550) 
ХРТ < ХРТ. + XINC 
GO TO 555 


DO 332 L=1,300 


J J*1 
R RET 
BCJ) = 
ACJ) = 


YPT 
(-CCXB( J)-DD)/BB 
IF (TABLE .NE. YES) GOTO 552 
WRITEC6,552) ACJ), BCJ), SIGMACM) 
IF (R410. - J/10) 332,331,332 
CALL ROOTS (ACJ), BCJ), AJ, BJ, CJ, N02) 
CALL EXCMS C'CLRSCRN') 
-WRITEC6,550) 
YPT = YPT + XINC 


CALL EXCMSC'CLRSCRN') 

IF (GRAPH .EQ. YES) CALL PLOTD(CA,B,J, .FALSE., 
LABEL, 'ALPHAS','BETAS', 
MINMAX,' S-$',SIGMACMD, 
XMIN, XMAX, YMIN, YMAX, GRD) 

IF (GRAPH .EQ. YES) GRD = GRD*1. 

NTINUE 


UE 


оз >> 0 <. 


Cn" Il 
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CONSTANT ZETA-OMEGA PLOTS 
IF (NZW) 359,359,350 


суо 


550 CALL EXCMSC'CLRSCRN') 
IF (TABLE .EQ. YES .OR. GRAPH .EQ. YES) WRITE (6,553) 
XWN = WN 


DO 358 M=1,NZW 


fe (TABLE .EQ. YES) WRITE (6,554) 
- 0 
R = 0. 
Я AZETA = 1.7300. 
DO 356 L=1,300 
ХИЧ = ZW(MO/AZETA 
D1 = 0.0 
D2 = 0.0 
С1 = 0.0 
С2 = 0.0 
В1 = 0.0 
В2 = 0.0 
C 
DO 353 N=1,NC 
K = N-1 
IF (K) 352,351,352 
351 91 = 0.0 
Q = -1.0/XMWHNxx2 
552 D2 = CJCN)¥Q1+D2 
C2 = BJCN)*Q1+C2 
B2 * AJCNOXQ1-*B2 
D1 = CJ(N)XQ+D1 
Cl = BJ(N)XQ+C1 
Bl = AJ(N)XQ+B1 
Q2 * -2.0XxZWCMDOXQ1-XWHNXx2xQ 
Q = 01 
ыы 91 = 92 
ТЕ (АВ5(В1*С2-В2%С1)-2) 356,356,354 
554 J = J*1 
R = К+]. 
ACJ) = C(C1xD2-C2xD1)Z(B1xC2-B2xC1) 
4 BCJ) = CB2XD1-B1XD2)/CBl1XC2-B2XC1) 


IF (TABLE .NE. YES) GO TO 356 
WRITEC6,552) ACJ), BCJ), ZWCM) 
IF (R710. = J710) 356,355,356 
555 CALL ROOTS (ACJ), BCJ), AJ, BJ, CJ, N02) 
CALL EXCMSC'CLRSCRN') 
WRITE (6,554) 
се AZETA = AZETA+(1.7300.) 


CALL EXCMSC'CLRSCRN') 


e 
IF €J .GT. 0) GOTO 357 
ИКТТЕС6 , 549) 
GOTO 558 
С 
557 IF (GRAPH .EQ. YES) CALL PLOTDCA,B,J,.FALSE., 
X LABEL, 'ALPHAS', 'BETAS', 
X MINMAX, 'ZW-$',ZWCMD, 
2 X XMIN,XMAX, YMIN, YMAX, GRD) 
IF (GRAPH .EQ. YES) GRD = GRD+1. 
e 
558 CONTINUE 
C 
7 CONTINUE 


100 


ΩΩ 


O00000 O O 


CONSTANT OMEGA PLOTS 
375 IF (NW) 385,385,376 
376 CALL EXCMSC'CLRSCRN') 
IF (TABLE .EQ. YES .OR. GRAPH .EQ. YES) WRITE (6,555) 


ро 584 М<1, МИ 


IF (TABLE .EQ. YES) WRITE (6,547) 
= Ü 
R = 0. 
AZETA = 0.0 
DO 382 L=1,300 
р1 = 0.0 
D2 = 0.0 
SL = 0.0 
С2 = 0.0 
В1 = 0.0 
B2 - 0.0 
DO 379 N=1,NC 
© K = N-1 
IF (K) 378,377,378 
377 U = 0.0 
Ul = -1.0 
378 U2 = 2.0XAZETAXU-U1 
D1 = (-1.00XXKXCJ (N) XHCM) XXKXUL-*DI 
D2 = (-1.00XxxKxCJ (N)XMHCM)XxxKxU*D2 
Cl = (-1.00XxKXBJ CN)XMHCMOXXKXUL4CI 
C2 = (-1.0)XXKXBJ (NDXHCM)XXKXU-*C2 
Bl - (-1.0)XXKXAJ CNO)XHCM)XXxKXUL*B1 
B2 * (-1.00XXKXAJ (NO) XMCM)XxxKXU- B2 
Ul = U 
Ой 2 
579 CONTINUE 
IF (ABS(CB1XC2-B2XC1)-2) 582,582,580 


380 J = Jtl 

R = Rtl. 

are = RARO AER ст: 
B 


B2XD1-B1xD2)/(BlxC2-B2xC1) 


IF (TABLE .NE. YES) GO TO 382 
WRITEC6,548) ACJ), BCJ), WCM), AZETA 
IF (R410. - J410) 382,381,382 

381 CALL ROOTS (ACJ), BCJ), AJ, В), СЛ, МО) 
CALL EXCMSC'CLRSCRN') 
WRITE (6,547) 

382 AZETA * AZETA*(C1./299.) 
CALL EXCMSC'CLRSCRN') 


IF (J .GT. 0) GOTO 383 


WRITEC6,549) 
GOTO 384 
585 IF (GRAPH .EQ. YES) CALL PLOTDCA,B,J,.FALSE., 
X LABEL, "А|РНАз", "ВЕТАз", 
Х MINMAX,' MW-$',MWCMD, 
X XMIN, XMAX, YMIN, YMAX, GRD) 
584 CONTINUE 
585 CONTINUE 
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ОО 


ОО 


ОО 


ОО 


400 


401 


402 


403 


404 


405 


>< >< >< 


IF (GRAPH „ЕФ. YES) CALL PLOTD(0.,0.,0.. .TRUE., 
f т Li т т T 
0.,! $',-9.7531, 
0.,0.,0.,0.,GRD) 
IF (GRAPH .EQ. YES) GRD = GRD+t1. 
GRD = 0. 
CALL EXCMSC'CLRSCRN'!) 


IF (OPT .EQ. YES) CALL DONEPL 

ТЕ СОРТ „ЕФ. YES) STOP 

МКІТЕ(6,557) 

МКІТЕ(6,556) 

МКІТЕ(6,557) 

МКІТЕ(6,558) 

МКІТЕ(6,559) 

WRITEC6,560) 

МКІТЕ(6,561) 

МКІТЕ(6,562) 

WRITEC6, 563) 

WRITEC6, 557) 

CALL READI (IANS) 
IF (IANS .GT. 6 .OR. IANS .LT. 1) GOTO 400 
GO TO (100, 404, 401, 402, 403, 405) TANS 


ROOT FINDER OPTION 
CALL EXCMSC'CLRSCRN') 
WRITEC6, 5649) 
WRITEC6, 565) 
CALL READR (ALPHA) 
HRITE(6,566) 
CALL READR (BETA) 
CALL ROOTS (ALPHA, BETA, AJ, BJ, CJ, N02) 
GO TO 400 


SAVE OPTION 
CALL SAVIT 
GO TO 400 


CREATE DISSPLA METAFILE OPTION 

CALL EXCMSC'CLRSCRN') 
МКІТЕ(6,567) 
МКІТЕ(6,568) 
МКІТЕ(6,569) 
CALL READC (ОРТ) 
IF (OPT .NE. YES) GOTO 400 

CALL. EXCMSC'CLRSCRN') 

WRITEC6,570) 

CALL DONEPL 


COMPRS = SUBROUTINE TO LET DONEPL FINISH 
CALL META 
GO TO 211 


SAME PROBLEM OPTION 
МКІТЕ(6,571) 
CALL READI CMINMAX) 
IF (MINMAX .EQ. 1) GO TO 200 
CALL EXCMSC'CLRSCRN'!) 
МКІТЕ(6,572) 
CALL READR (ХМІМ) 
WRITEC6,573) 
CALL READR (XMAX) 
ИКТТЕ(6, 574) 
CALL READR CYMIN) 
WRITEC6,575) 
CALL READR CYMAX) 
GO TO 200 
CONTINUE 
RETURN 
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500 


501 
502 
505 
504 
505 
506 
507 
508 
509 
510 


511 
512 
513 
514 
515 


516 


518 
519 
520 
521 
522 
525 
524 
525 
526 
527 
528 
529 
550 
551 
552 


534 
535 
536 
537 
538 


540 
541 
542 


FORMATS STATEMENTS 
FORMATC5X,'THIS IS THE INTERACTIVE PARAMETER PLANE PROGRAM...',/, 
5X,'THE USER WILL BE PROMPTED FOR VARIOUS INPUTS.',///7, 
5X,'WILL YOU BE ENTERING DATA FROM A CONSOLE OR DATAFILE?', 
/Z,15X; етут ОР WC) 
FORMAT(//^/,' ENTER TITLE TO APPEAR FOR THIS FAMILY OF CURVES') 
FORMATC/,1X,'WHAT IS THE STARTING VALUE OF OMEGAN CWN>0.0)7') 
FORMAT(/,1X,'WN MUST BE GREATER THAN ZERO - TRY AGAIN') 
FORMATC/,1X,'HOM MANY DECADES PAST WN ARE DESIRED? (ND)') 
FORMAT(C/,' МИНАТ IS THE ORDER OF THE CHARACTERISTIC EQUATION?(CNO2') 
FORMAT(C/,1X,'HOM MANY CONSTANT ZETA CURVES ARE DESIRED? (NZ)') 


FORMAT(/,' ENTER THE VALUES OF ZETA TO BE USED IN COMPUTATION. ..') 
FORMATC5X, 'ZETAC',I2,') = 7°) 
INCLUSIVE - TRY AGAIN") 


FORMATC5X,'ZETA MUST LIE BETWEEN 0 AND 1, 
HOW MANY CONSTANT SIGMA CREAL ROOT) CURVES ARE DESIRED? 


ҒОКМАТ(/,! 
X (NS)') 
FORMAT(5X,"SIGMA(',I2,") г 
FORMATC5X,'NEGATIVE SIGMA MEANS POSITIVE REAL ROOT - TRY ANOTHER') 
FORMATC/,1X,'HOW MANY CONSTANT WN CURVES ARE DESIRED? СМИ) ') 
FORMATC5X, 'WC',I2, ') 24) 
ҒОКМАТ (5X, TWN NOT WITHIN PLOTTABLE RANGE’, 
/ 5X, "YOUR USABLE RANGE До 
/,10Х,Е10.2,' TO ',F10.2) 
FORMAT(/, 1X,"HOW MANY CONSTANT ΖΧΥΙΝ CURVES ARE DESIRED? (NZMW2') 
FORMATC5X, 'ZNC',I2,') < 7%) 
FORMAT(5X,'NON-POSITIVE Z-WN MEANS POSITIVE ROOT - TRY ANOTHER') 
FORMATC/,1X,'ENTER THE CONSTANT COEFFICIENTS...') 
FORMAT(C5X,'---Sxx',I2,' CJC', I2, '0 = 7°) 
FORMATC/,1X,'ENTER THE ALPHA COEFFICIENTS...') 


X 
X 
X 


X 
X 


EBRMAT(5X,'---5XX',I2,'--- А.,/(!,12,") - 7%) 

ҒОКМАТ(/,1Х, "ЕМТЕК THE BETA COEFFICIENTS...') 
FORMAT(5X,'---SXX',I2,'--- BJC',I2,") = 7°) 

FORMATC/,' WANT TO REVIEM YOUR ENTRIES BEFORE RUNNING? (Y^/N2') 
FORMAT (/,10X,' GRAPH TITLE') 

FORMAT (1X,9A4) 

FORMAT (7,8X,2HND, 8X, 2HNO, 8X, 2HNZ, 8X, 2HNS, 8X, 2HNW, 7X, SHNZW) 
FORMAT (6110) 

FORMAT (7,10X,"INITIAL VALUE OF OMEGA = ',F10.5) 

FORMAT (/,10X,' ZETA ') 

FORMAT (8E10.3) 

FORMAT (1X,'..... NO VALUE..... b) 

FORMAT (/,10X,'SIOMA ὁ) 

FORMAT (/,10X,' W D) 

FORMAT (/,10X, 'ZM Те) 

FORMAT (/,10X,'CONSTANT COEFFICIENTS IN DECENDING ORDER') 
FORMAT (/,10X, 'ALPHA COEFFICIENTS IN DECENDING ORDER') 
FORMAT (/,10X,'BETA COEFFICIENTS IN DECENDING ORDER') 
FORMAT (/,10X, 'XMIN XMAX YMIN ҮМАХ") 

FORMAT (1X,4E10.3) 


FORMATC/,' WANT TO MAKE ANY CHANGES? CY/N)') 


ҒОКМАТ(/, 
X5X, ' TITLE 


" WHAT VARIABLE“AREA DO YOU WISH TO CHANGE?',//, 
TI OMEGA START..WN * DECADES....ND',/ 
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566 
5@5 
546 


548 
549 


550 


552 
555 
556 


556 
557 
558 


560 
561 
562 
565 
566 
565 


566 
567 


568 


569 
570 


571 


57e 
57 3 


575 
576 


XBX, Όκρε ΟΡ М0 CONST ZETA...NZ CONST SIGMA..NS',/ 
X5X, "CONST WN..... NW CONST Z-MWN...ZW',77 

X5X, "CONST TERMS. .CJ ALPHA TERMS..AJ BETA TERMS...BJ',7/ 
ХХ ЕНИ САТТЕ Е NEN PROBLEM..PR NO CHANGE....NC',///, 
X1X, "ENTER TWO DIGIT CODE...') 


FORMATC/,' DO YOU WANT TABULATED DATA ON THE SCREEN? (Y/N)') 
FORMATC/,' DO YOU WANT THE GRAPHS ON THE TERMINAL? CY/N2') 
FORMAT C1H1,10X,'CONSTANT ZETA CURVES'!) 
FORMAT(/,10X, 'ALPHA BETA OMEGA ZETA') 
FORMAT (€4E16.5) 

ΧΤΗΝ аи DUE TO PLOT RESTRICTIONS, A COMPLETE GRAPH CANNOT BE OU 
FORMAT C1H1,10X,* CONSTANT SIGMA CURVES") 
FORMAT (7,10X, "ALPHA BETA SIGMA') 
FORMAT (3E16.5) 
FORMAT C1H1,10X,'CONSTANT ZETA-OMEGA CURVES') 


FORMAT (/,10X, "ΑΙ ΡΗΑ BETA ZETA-OMEGA ') 
FORMAT (1H1,10X,'CONSTANT OMEGA CURVES') 
FORMATC5X,'| OPTION NO.| OPTION |!) 
ҒОКМАИТПСБА, !----------------------2.......4.. 22-22 0050 |") 
ҒОКМАТ(5Х, ! | 1 | NEW PROBLEM Ы. 
FORMAT C5X, | 2 | SAME PROBLEM |!) 
FORMAT C5X, | 3 | ROOT FINDER |") 
FORMAT C 5X, | 4 | SAVE DATA |!) 
FORMAT(5X, ' | 5 | SAVE GRAPH IN DISSPLA METAFILE |!) 
FORMATC5X, ' | 6 | RETURN TO MAIN MENU |") 
FORMATC5X, 'ENTER AN ALPHA-BETA PAIR, AND THE ROOTS OF YOUR τ; 
Х 5X,'SYSTEMS CHARACTERISTIC EQUATION HILL BE RETURNED г) 
ҒОКМАТ(///,5Х, "ENTER THE ALPHA VALUE та 
ҒОКМАТ(/,5Х, "ЕМТЕК ТНЕ ВЕТА VALUE Ὁ) 


FORMATC5X,'YOU NOW HAVE THE OPTION OF STORING THE LAST SET OF  '7, 


X 5X,'CURVES IN A DISSPLA METAFILE. THIS ALLOWS RETRIEVAL'/, 
X 5X, 'OF DATA AT A LATER TIME FOR ROUTING TO ANY OF "Ж 
X 5X,'SEVERAL OUTPUT DEVICES CTEK618, 3800 LASER PRINTER, "и, 
X 5X,'VERSATEC PLOTTER, ETC.) 1//) 

FORMATC5X,'IF YOU CHOOSE THIS OPTION, THE PROGRAM MUST BE K 
X 5X,'TERMINATED - THIS CANNOT BE AVOIDED WITHOUT n 
X 5X,' 'CATASTROPHIC RESULTS. '// 79 

FORMAT(///,5X,'DO YOU WISH TO USE THIS OPTION? τ 
Χ 154; "ΚΝ. τν 72 

FORMAT(/////,5X,'IF YOU WISH GRAPHIC OUTPUT, TYPE: 7, 
X 15X,' '"DISSPOP" '/ 
X 5X,*AND FOLLOW THE INSTRUCTIONS . . E 


FORMAT (//7,' AUTOSCALE OR USER-DEFINED LIMITS FOR CURVES?', 
Хи," 1L=AUTOSCALE; O=USER-DEFINED') 

FORMAT (//,' INPUT MINIMUM VALUE FOR X (X-MIN2'!) 

FORMAT (/,' INPUT MAXIMUM VALUE FOR X (Х-МАХ)") 

FORMAT Си," INPUT MINIMUM VALUE FOR Y (СҮ-МІМ)") 

FORMAT (/,' INPUT MAXIMUM VALUE FOR Y CY-MAX)') 

ΕΟΚΜΑΤ (η) 

ΕΝΡ 
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SUBROUTINE PLOTD -- GRAPHS WILL BE PRODUCED ON UPRIGHT 11 X 14 
PAGE WITH 9 INCH AXES. IF USER SELECTS 'AUTOSCALE' FEATURE, 
SUBROUTINE PLD010 CINTERNAL TO PLOTD) FINDS MIN AND MAX FOR EACH 
AXIS AND SCALES ACCORDINGLY. FORMAT: 


CALL PLOTDCXDATA, YDATA, NNPTS, EJECT, LABEL, XLABEL, YLABEL, 
X 


WHERE: 
XDATA 


YDATA 


NNPTS 


EJECT 


LABEL 


XLABEL 


YLABEL 


MINMAX 


CRVTTL 


CRVNUM 


MINMAX, CRVTTL, CRVNUM, XMIN, XMAX, YMIN, YMAX) 


IS A REALX4 ARRAY DIMENSIONED AT LEAST |NNPTS| 
CONTAINING THE X ORDINATE VALUES, 


IS A REAL*¥4 ARRAY DIMENSIONED AT LEAST | ММРТ5 | 
CONTAINING THE Y ORDINATE VALUES, 


15 ΑΝ INTEGER¥4 SCALAR DESIGNATING THE NUMBER OF 
POINTS TO BE PLOTTED. THE NUMBER OF POINTS IS 
ABSCNNPTS). NNPTS«0 MEANS PLOT POINTS ONLY. 


IS A LOGICALX4 VARIABLE OR CONSTANT INDICATING 
WHETHER A PAGE EJECT IS REQUIRED FOLLOWING THE 
CURRENT CURVE. THIS ALLOWS MULTIPLE CURVES ON ONE 
SET OF EXES. PAGE EJECT WILL OCCUR FOR NEXT GRAPH 
AFTER EJECT HAS BEEN SET TO .TRUE. 


IS A QUOTED LITERAL OR HOLLERITH STRING OR ARRAY 
CONTAINING THE INTENDED LABEL FOR THE GRAPH. THE 
MAXIMUM ALLOWABLE LENGTH CINCLUDING '$' CHARACTER) 
IS 32 CHARACTERS. 


IS A QUOTED LITERAL OR HOLLERITH STRING OR ARRAY 
CONTAINING THE INTENDED LABEL OF THE X-AXIS OF THE 
GRAPH. IN THIS PROGRAM, XLABEL IS ALWAYS ‘ALPHA’. 


IS A QUOTED LITERAL OR HOLLERITH STRING OR ARRAY 
CONTAINING THE INTENDED LABEL OF THE Y-AXIS OF THE 
GRAPH. IN THIS PROGRAM, YLABEL IS ALWAYS 'BETA'. 


IS A PARAMETER THAT DETERMINES WHETHER THE MINIMUM 
AND MAXIMUM VALUES FOR THE AXES ARE TO BE ASSIGNED 
BY THE USER, OR WHETHER THEY WILL BE ‘AUTOSCALED'. 


IS A QUOTED LITERAL OR HOLLERITH STRING OR ARRAY 
AND TERMINATED BY А 9! CHARACTER SPECIFYING THE 
INTENDED NAME WHICH LABELS AN INDIVIDUAL CURVE. 


IS A REAL VARIABLE OR CONSTANT THAT SPECIFIES THE 
VALUE TO BE CONCATENATED ONTO THE END OF '"CRVTTL'. 
FOR EXAMPLE, IF THIS CURVE REPRESENTS ‘ZETA = 0.5' 
THEN CRVTTL = 'Z=$', WHILE CRVNUM = 0.5. 


w < < < αν κα s < s αν αν ω ϱ-- αν αν ee ie eee ee ie s ee ieee ei i ee ee ee αν καν αν αν αν αν αν καν καν αν καν αν αν αν καν καν καν «мәз «шо еш» 
- - HH Hsu JH αν καν чай чо απ 


SUBROUTINE PLOTDCXDATA, YDATA, NNPTS, EJECT, LABEL, XLABEL, 
X 


X 
КЕА| ха 
REAL xà 
REAL ¥4 
INTEGERX4 
LOGICAL *4 
LOGICALXI1 
LOGICALXI1 
LOGICALXI1 


( ОбІСАІ Χά 


YLABEL, MINMAX, CRVTTL, CRVNUM, 
XMIN, XMAX, YMIN, YMAX,GRD) 
XDATAC1), YDATA(C1) 
XMIN, XMAX, YMIN, YMAX 
CRVTTL,CRVNUM 
NPTS,NNPTS 
EJECT 
LABELC1) 
ΧΙ ΑΒΕΙ (1) 
ΥΙΑΡΕΙ (1) 


INIT /.FALSE.^/ 
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O 


о О 


10 


20 


50 
40 


~ Сх 
оо 


80 


SET THE NUMBER OF POINTS 
NPTS = IABSCNNPTS) 


IF THE ROUTINE HAS BEEN NOT BEEN INITIALIZED CINIT = .FALSE.) 


THEN INITIALIZE IT. 


IF €.NOT. INIT) CALL PLDOOICXDATA, YDATA, NPTS, LABEL, XLABEL, 
Х 


YLABEL, MINMAX, CRVTTL,CRVNUM, 
XMIN, XMAX, YMIN, YMAX,GRD) 


FRAME THE PLOT AND REORDER SYMBOLS 
IF €.NOT. INIT) CALL FRAME 


INDICATE INITIALIZATION IN CASE THIS ROUTINE IS RE-ENTERED 
WITH MULTIPLE CURVES. 


INIT = .TRUE. 


E за THE NUMBER OF POINTS TO GIVE ABOUT 5 MARKERS PER 
LINE 


TENTATIVELY SET POINTS ONLY, THEN CHECK 
NMARK = -1 
IFCNNPTS.LE.02GO TO 10 


CURVE WANTED, SET MARKERS 

IF (MODCNPTS, 4) .EQ. 1) NMARK 
IF CMODCNPTS, 4) „МЕ. 1) ММАКК 
IF CNMARK .EQ. 0) NMARK 7 1 
CONTINUE 


NPTS / 6 
МРТ5 / š 


DRAW THE CURVE 


XXMAX -1.0Е75 
ҮҮМАХ =1.U0E75 


DO 20 I-1,NPTS 
IFCXDATACI) .GT. XMAX .OR. XDATACI) .LT. XMIN „ОК. 
YDATACI) .GT. YMAX .OR. YDATACI) .LT. YMIN) GO TO 20 
IFCXXMAX .LT. XDATACI)) J-I 
IiFCOXXMAX .LT. XDATACIO) XXMAX < ХОАТА(Т) 
IFCYYMAX .LT. YDATACI)) K-I 
IFCYYMAX .LT. YDATACI)) YYMAX = YDATACT) 
CONTINUE 


ІҒССҮМАХ-ҮҮМАХ) - (XMAX-XXMAX)) 40,30,30 

L= 

60 ТО 50 

L=K 
CALL GRACECO.2) 
CALL DOT 
IF (GRD) 60,60,70 
CALL GRID(1,1) 
CALL RESETC'DASH') 

IF CCRVNUM „Ед. -9.7531) GO TO 80 

CALL HEIGHTCO.125) 
CALL RLMESSCCRVTTL, 3, XDATACLO, YDATACLO)D 
CALL RLREALCCRVNUM, 2, ' ABUT' , ' ABUT !D 
CALL THKCRVC0.015) 


CALL CURVECXDATA, YDATA, NPTS, NMARK) 
CALL RESETC'THKCRV'!) 


IF THIS IS NOT THE LAST COR ONLY) CURVE ON THIS GRAPH, THEN 
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EXIT. OTHERWISE CLOSE THE PLOT AND TURN OFF INITIALIZATION 
FLAG. 


IF C. NOT. EJECTO RETURN 
END OF THIS PLOT 

CALL ENDPL(O) 

INIT = .FALSE. 


RETURN 
END 


SUBROUTINE PLDOOI1CXDATA, YDATA, NPTS, LABEL, XLABEL, 
X 


YLABEL, MINMAX, CRVTTL, CRVNUM, 
XMIN, XMAX, ҮМІМ, YMAX, GRD) 


THIS SUBROUTINE DOES THE INITIALIZING. 


REAL X4 XDATACNPTS) 
REAL xà YDATACNPTS) 
INTEGERX4 NPTS 
LOGICALX1 LABEL(1) 
LOGICAL *1 XLABEL(1) 
LOGICAL *1 YLABELC1) 


ΚΕΑΙ Χά ХМІМ 
КЕЛДІ +4 ХМАХ 
ΚΕΑΙ Χά ΥΜΙΝ 
REAL Xà YMAX 


INITIALIZE DISSPLA 


CALL PLD009 

CALL HEADINCLABEL, 100, 
CALL XNAMECXLABEL, 100) 
CALL YNAMECYLABEL, 100) 


2, 1.) 


EXTRACT MINIMA AND MAXIMA 

IF (МІММАХ .МЕ. 12 60 Т0 90 

CALL PLDOIOCXDATA, МРТ5, ХМІМ, ХМАХ) 

CALL PLDOLOCYDATA, NPTS, YMIN, УМАХ) 
CALL THE LINEAR-LINEAR INITIALIZING ROUTINE 

CALL PLDO11 CXMIN, XMAX,YMIN, YMAX) 

RETURN 

END 

SUBROUTINE Ριρ009 
THIS SUBROUTINE ESTABLISHES THE PARAMETERS FOR DISSPLA. 
NOTE THAT IT IS THE USER'S RESPONSIBILITY TO NOMINATE THE 
GRAPHIC DEVICE. 


CALL NOCHEK 

CALL NOBRDR 

CALL PAGE(14.,14.) 

CALL PHYSOR(2.,.75) 
GO TO LEVEL 2. 


CALL AREA2D(9.,9.) 
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LETTERING IS DUPLX WITH UPPER CASE ONLY. 


CALL DUPLX 
CALL BASALFC'STANDARD') 


INTEGER COR ROUNDED) AXES WITH Y AXIS LABELLING AT 0 DEGREES. 


CALL INTAXS 
CALL YAXANGCO.) 
RETURN 

END 


ΩΩ 


ООО 


SUBROUTINE PLDO1OCV, N, MIN, MAX) 
THIS SUBROUTINE SCANS A VECTOR FOR MAXIMUM AND MINIMUM 
INPUT PARAMETERS: 


V DATA VECTOR (REAL) 
N NUMBER OF POINTS IN VECTOR V (INTEGER) 


OUTPUT PARAMETERS: 


MIN VECTOR ORIGIN (REAL) 
MAX VECTOR MAXIMUM (REAL) 


REAL X4 VCND 
INTEGERX4 N 

REAL xà MIN 
REAL x4 MAX 


INITIALIZE THE MAXIMA AND MINIMA 


I SUET5 
-1.0E75 


FIND MAXIMUM AND MINIMUM OF VECTOR V 


DO 100 I = 1, N 

IF (MIN .GT. VCI)) MIN 

IF (MAX .LT. VCI)) MAX 
100 CONTINUE 


(3 (3 (> (3) С) С) (3 С) С? (> (3 (КФ) (300 


ООО 


МІМ 
МАХ 


ООО 
nu 


поп 
<< 


(1) 
(1) 


RETURN 
END 
Е 
С SUBROUTINE PLDO11 CXMIN, XMAX, YMIN, YMAX) 
с THIS SUBROUTINE SETS UP DISSPLA FOR A LINEAR-LINEAR AXIS PLOT. 
REAL xà XMIN 
REAL x4 XMAX 
REAL X4 YMIN 
4 ΚΕΑΙ Χά ΥΜΑΧ 
с A SIMPLE CALL TO GRAF WILL DO IT... 
CALL HEIGHTC0.175) 
CALL GRAFCXMIN, "SCALE", XMAX, YMIN, "SCALE", YMAX) 
RETURN 
END 
C 
С 
С 
С 
С 
с 
C 
С 
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SUBROUTINE ROOTS -- CALCULATES ROOTS OF THE NO ORDER EQUATION 


FOR EVERY TENTH VALUE ALPHA/BETA PAIR. 
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SUBROUTINE ROOTS CALPHA, BETA, 
INTEGER¥4 NO2, NN, NNN, NNNN 

REALX6 ALPHA, BETA, AJ(100), 

REALX8 COEF(C(10050, ROOTMY(100) 
WRITEC6,40) 

NN = N02+1 

NNNN = NO2 

DO 10 I-1,NN 

NNN - NN*1-I 


AJ, B 
BJC100) 


10 COEFCNNN) = CJCI) + CALPHA x AJ(I)) + (BETA x BJ(I)) 


CALL ZPOLR (COEF,NO2,ROOTMY,IER) 
WRITE (6,50) 
DO 20 1<1, ММММ 
II - 2XI 
M= II-l 
20 WRITE (6,60) ROOTMY(III), ROOTMYCIT) 
WRITE (6,30) 
30 FORMAT(//7/4/444444444444 447) 
40 ҒОКМАТ(/20Х," ROOTS FOR ABOVE ALPHA, 
50  FORMATC21X,' REAL PART 
60 FORMATC21X,E10.4,6X, E10.4) 
RETURN 
END 


J, 


CJ, 


N02) 
, CJ(100) 


BETA') 


IMAG PART') 


SUBROUTINE SAVIT -- SAVES DATA IN FN FT 
WHERE INAME IS THE USER'S CHOICE. 
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SUBROUTINE SAVIT 
COMMON ZSAVEZ LABEL, WN, ND, 
NZ, ZETA, NS, 
XMIN, XMAX, 
КЕАС жа NN, CJ(100), 
КЕЛІ жа ФЕТАС100), УТСМАС100), ИМ(100) 
INTEGER ND, NO2, NC, NZ, 
CHARACTERXG4 LABEL(9), 
МКІТЕ(6,10) 
READC5,20) CINAME(CID, 


мог, NC 


>< >< 


YMIN, YMAX 
AJ(100), BJ(100) 


ІМАМЕС2) 
I-1,2) 


NS, NW, NZW 


CALL FRICMS 
WRITEC2,30) 
WRITEC2, ®) 
WRITEC2, x) 
WRITEC2,X) 
WRITEC2,X) 
WRITEC2,X) 
WRITEC2,X) 
WRITEC2,X) 
ИКТТЕС2,%) 
ИКТТЕС2,%) 
WRITEC2,X) 


10 FORMATC5X,'UNDER WHAT NAME DO YOU WANT TO SAVE THE DATA? ',7/, 
X БХ, "(8 CHARACTERS MAX) ') 


20 ҒОКМАТ(2А4) 
50  FORMAT(C9A4) 
END FILE 02 
REWIND 02 
RETURN 
END 


C'FILEDEF ','02 
CLABELCI), I=l, 
WN 
ND, NO2, NC, NZ, 
(CJ(J), J=1, NC) 
(AJ(J), J=1, NC) 
(BJ(J), J=1, NC) 
(ZETA(M), M=1, 
CSIGMACM), ΜΞΙ, 
(WCM), M=1, NW) 
(2И(М), M=1, 
ХМІМ, ХМАХ, 
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9) 
NS, 


NZ) 
NS) 


NZW) 
YMIN, 


vU T) 


NW, 


YMAX 


» С 


SIGMA, NW, 
» XMIN, 


J; 
W, 


AJ, 
ΝΖΗ, 


, ΖΗ(100) 


І5К 


NZW 


", INAME, 'DATA 


BJ, 
ZW, 


XMAX, 


YMIN, YMAX 


00000000000 


 — æ Åe «әлә «ӘӘ «ӘӘ ә ӘӘ s“ s“ G s“ s“ G G s G “G G s G s 4“ G G G G s“ G G Q s s < Q Q G s s шый = s «жә m 2 «ни что «то «ню «е «но = «αν | «πα «πὶ 
ж: «на «но сиз «к. «ж. «ж. «ж. «ко s «ню «жә сию ο «ню s αν αν αν ο s“ «но s «ню s s G“ s s s G G“ G G «ню G G G < G G s G s G G G G ЧЫ G G ШЕР шый ыў ШЕ ыз s G s s G s s s «πο s s s 


SUBROUTINE GETIT -- RETRIEVES DATA FROM FN FT FM = INAME DATA Al, = 
WHERE INAME IS THE USER'S CHOICE. = 


t Q “b Qm тиў чый шыў «и «αὶ ο am чый} чыў ЧЫ ыў ay GP aw aw | « «жә «4 « G G « G G G тый сз сз ср “ик: «жә «но «но «ню «но G «кз «ы; πο απο «жә «ы ему «εἴν G s G G G G“ «не «жә «жа s «не «не «жә απο απο «но «но s“ 
ж. «нә чый «ы сый шый ϱ «ы «жә «то «ж: «жә «но «ний жыў ϱ  Ὢ ο ο «Ὢ «α «нә «жз «но «из G G G G G G G G G «но «но»  - - -. -- «из «9 тыў шый тый «не s“ s. «но s «но s. s“ s «на «на «не «из «иә «не «нә «нә «не «шә 


10 


20 


50 
40 


50 
60 


>< >< 


SUBROUTINE GETIT 
COMMON /SAVE/ LABEL, WN, ND, NO2, NC, CJ, AJ, BJ, 
NZ, ZETA, NS, SIGMA, NW, W, NZW, ZW, 
ХМІМ, ХМАХ, ҮМІМ, ҮМАХ 
КЕА1*6 ИМ, С./(100), А.(100), В./(100), ХМІМ, ХМАХ, ҮМІМ, ҮМАХ 
КЕЛІ ж4 ZETAC100), SIGMAC100), WC100), ZWC100) 
INTEGER ND, NO2, NC, NZ, NS, NW, NZW 
CHARACTERX4 LABEL(C9) 
СНАКАСТЕКЖ8 ІМАМЕ 
CHARACTERX21 NAME 
МКІТЕ(6,40) 
READC5,50) INAME 
NAME = "STATE '//INAME//' DATA x' 
CALL EXCMSCNAME, RC) 
IF (RC .EQ. 0) GOTO 20 
WRITEC6, 30) 
GOTO 10 
CALL FRTCMSC'FILEDEF ','02 ", "DISK ", INAME, "DATA г) 
READC2,60) CLABELCIO, I-1, 9) 
READC2,X) ИМ 
READC2,X) ND, N02, NC, NZ, NS, МИ, МЕЙ 
READC2,X) (CJ(J), J=1, NC) 
READC2,X) (AJ(J), J=1, NC) 
READ(2,%x) (BJ(J), J=1, NC) 
READ(2,X) (ZETA(M), M=1, NZ) 
READ(2,%x) (SIGMA(M), M=1, NS) 
READC2,X) (ИСМ), М=1, МИ) 
READC2,X) CZWCMD, M=1, NZH) 
READC2,X) XMIN, XMAX, YMIN, YMAX 
FORMATC5X,'DATA FILE NOT FOUND - TRY ANOTHER') 
FORMATC5X, 'UNDER WHAT NAME IS THE DATAFILE SAVED? ",/, 
5Х, "(8 CHARACTERS MAX)') 
FORMATC1A8) 
ҒОКМАТ(9А4) 
END FILE 02 
REWIND 02 
RETURN 
END 


«ж: «жо тый ЧЫР #0 «шм ы) тый шыў =} «то «тв G G G G G «нә «на «кв | « ЧЫР ыр тый чыў ый απο G G G G G «но G G G G G G απο «то «не «ж: «из «ко «кә «ка «ж: «ню guy rx 
— “< «и «4 «4 « «жә «жә «не «нә «ко: «ше ϱ сый еш, «кә s G G G s G G G G G G G G «но» G G G s G Qm Q G s G s G G G s απο «нә «но «но «но» «не «и «нә «нә «на «нә «из «шә 


SUBROUTINE ASTER(LLINES,LINE) 
CHARACTERX4 DOLLAR/'$ '/, LLINESC8)D, LINEC9) 
DO 10 I-1,8 
LINECID-LLINESCID 
CONTINUE 
LINEC9)=DOLLAR 
RETURN 
END 
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SUBROUTINE READR -- INTERACTIVELY READS A REAL NUMBER REPLY. - 
IF THE USER INADVERTENTLY ENTERS A NULL STRING i - 
A WARNING IS ISSUED AND ONE RECOVERY IS ALLOWED. = 


- - «нә «к « eee өне ee ==} s“ өне «нә «не ене «нә сир «πὶ s Gum «и eee «и! то «но что ee aD aa eee «ке «то (ШЫР md «и «πὶ ЧЕР «ы «и «καν «πὶ απο απο ты} «ын ты} сы mdp сы} сы} «и «и шыр s G шыр «и? ты} шр сто «но a= өше 
сно сир ap аар аар ано ано сар Әә «ы% s «и αν өне тыз тыз тыз өне «не тыз тыр «не απο απο απο αν αν Ер шыр cnp G G «и? өне өне «на «но «и G но Ер «нә «| « ан «нә «но что ат ен: ене «но «но «нә «нә «нә «нә «мә «ы απο «| απ «шо «нә се өн. но um 


0000000000 


SUBROUTINE READR CANSR) 
REALX4 ANSR 
INTEGER COUNT 
COUNT =0 
10 CONTINUE 
СОЦНТЕСОЦНТ+ 1 
IF (COUNT.LT.3) GO TO 20 
МКІТЕ (5,60): 
60 Т0 40 
20 CONTINUE 
READ (5,%,END=30,ERR=30) ANSR 
RETURN 
50 REWIND 5 
WRITE (5,50) 
GO TO 10 
40 CONTINUE 
STOP 
50 FORMAT CIX,' WARNING: NULL STRINGS ARE NOT ALLOWED, ENTER A NUMER 
XICAL VALUE. ') 
60 Ui C///,5X,* PROGRAM TERMINATION - TWO NULL STRINGS ENTERED! ') 
N 


ашо тыз ене ошо еше еше сыр ашр ο ο αν «ων ο αν ο αν «ο απο απο «и? Ер өне ене ене «нә «но ШЫР ЧЫР ене ене екы но πω «но чо ШЫР өне ене ене ШЕР ене πω ШЫР ене (БЕР «ке «шо «ко «же UND «шо «ке «кә «кә «ко «жә am am «но «но «к» «шә «шо «шә «шә a= = 
Чыр тыр ==} ы: ешю «и «жо «к» «ко өш: «то Å ШЫР «но өне «нә «но ене ана ене ешю енә ене ене но ене «но «нә ешю ешю ешю «ню «но «не сш ешю «нс «нә «кә өн: «ко «ко тыз чыз ан: чыз «то «то «жо ωρα 


SUBROUTINE READI -- INTERACTIVELY READS AN INTEGER REPLY. z 
IF THE USER INADVERTENTLY ENTERS A NULL STRING OR NEGATIVE VALUE = 
A WARNING IS ISSUED AND ONE RECOVERY IS ALLOWED. = 


ә «жо ек. ЧЫР өне ене ешю ешо «нә ШЫР «не ШЫР «нә ене ешю ешю «нә ешю ешо ене ешю ешю «нә ене өне өне «кә ешю ешю ешю ешю ене өше ШЫР «ню сыр ЧЫР «өне «но «ню «но ЕР жы} ты} сы} αν «κα сәр ашр сәр сәр «и s өне ШЕР сы} ЧЕР ЧЕР ee ==} өне еше G s mb m өше «αν 
ғ. «но «но ене «но «но «но ене ене «но «но «но απο «ο «ο αν «но «но απο απο сыр UND сыр - сир ты} cn ты} сы} чы} ты} чы} шыр «и чы} «ын сы) сты} G s s G ωρα 


0000000500 


SUBROUTINE READI (IANS) 
INTEGER COUNT,IANS 
COUNT=0 

10 CONTINUE 
COUNT =COUNT+1 
IF CCOUNT.LT.3) GO TO 20 
WRITE (5,70) 
GO TO 50 

20 CONTINUE 
READ (€5,%®,END=40,ERR=40) IANS 
IF (IANS) 40,30,30 

30 CONTINUE 
RETURN 

40 REWIND 5 
WRITE (5,60) 


GO TO 10 
50 CONTINUE 
STOP 
60 ΝΑΙ (1X,' WARNING: IMPROPER DATA ENTRY! ENTER A POSITIVE INTEG 
70 RAT (ZZZ,5X," PROGRAM TERMINATION - 2 IMPROPER DATA ENTRIES!') 
С 
С 
С 
С 
5 
ς 
с 
с 
с 
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ООООО ООО 


10 


20 


ООО ООО 


-— -- < s чыр s“ «жз «ше «ж: «жә G 4 G «жә «нә «кв «ше «не «жә  α «ЕР ew wee wee ew «но ND dab αν αν αν αν αν αν «κ ο αν αν αν αν «ο чыз ЧЫР ЧЫР ЧЫР сы} cm  --ωω amb | «ο «πὶ 
ено СЕР СЕР ФЕР чыз «кә «ж» ШЫР SB Bess es -.  , «ЕР SB ou SB ww еже ош» «к» «шо «кә «ж» «к» «кә ЧЫР ЧЫР ЧЫР ЧЫР G G G G G G SFT FT SFT FT G «жы G s G «ο απο G «жә «жә G s. s s“ «жә = = 


SUBROUTINE READC -- INTERACTIVELY READS A CHAR STRING REPLY. = 
('*YES' OR 'NO'). IF THE USER INADVERTENTLY ENTERS A NULL STRING z 
A WARNING IS ISSUED AND ONE RECOVERY IS ALLOWED. I = 


q < s < G «кә «ж» «шә «нә «шә G «шә “G «ж» «иә «шә G G G G “G G G G G G «ше «шә «шә ЧЫР «шә «и» «ыз «| αν FT «и» «и» «но «но s G s “G G V8 G G «не» Ss “s G G «но «жә G G G G G s“ «и? = = = 
aga ешо «жә s «шо «жә «жә s“ «жә ЧЫР ЧЫР ЧЫР s G ЧЫР “G G G G G ЧЫР ЧЫР «но «жә ЧЫР ЧЫР ЧЫР ЧЫР ЧЕЗ ЧЫР «πὸ Sie «и s G G G Q x < G G G G G  .ϱϱ -α - Әәә ϱ s ϱ.ῥ- 


SUBROUTINE READC (CANS) 
INTEGER COUNT 

CHARACTERX4 CANS 

COUNT + 0 

CONTINUE 

COUNT -COUNT*1 

IF CCOUNT.LT.5) GO TO 20 

WRITE (5,60) 

GO'TO 40 

CONTINUE 

REWIND 5 

READ (€5,70,END=30,ERR=30) CANS 
RETURN 

REWIND 5 

WRITE (5,50) 

GO TO 10 

CONÍINUE 

STOP 

FORMAT C1X,' WARNING: NULL STRINGS ARE NOT ALLOWED ') 
FORMAT (///,5X,' PROGRAM TERMINATION - TWO NULL STRINGS ENTERED! 
Coe (A2) 


- en a a ФЕР сиф ФЕР ФЕР сәф «ο αν «πὸ απὸ αν απὸ απ" «πὸ «πὸ απ απὸ απὸ G“ сыф G G G G «ию «жә «но» απ ЧЫР G ЧЫР ЧЫР чыз «жә “s G G G «жә «ж. «к: G G G «ж» «жә G ЧЫР G G G ЧЫР G «ко «жә G s“ «ко «ез «шә «шә 
«кз «жс «шо «жо ЧЫР ЧЫР ЕР Nl  ϱ Be se G G «шо «ж «ж» «шә «шә «нә «жә «ж» «и» G G G G СЫР G G чие G «ы» чи: G G G G G G G G чыр «ик «и» «ж: «жә «жә «ыл» «с> «шо «жә «шә «шә 


SUBROUTINE READL -- INTERACTIVELY READS A STRING OF CHARACTERS. = 
IF THE USER INADVERTENTLY ENTERS A NULL STRING = 
A WARNING IS ISSUED AND ONE RECOVERY IS ALLOWED. = 


— a «πὸ 
ено ешо «шо 


10 
20 


50 


40 


50 
60 


70 
80 


=a Q «и? ο αρ «ο απὸ απὸ απ απ «ο s“ G «и» ШЫР БЫР See ee ЧЫР ЧЫР ЧЫР m E «шә «ез «жә ЧЫР G “G G «нә «нә «шә «ав αν «ο «| «| αν «ο «πὶ απ αρ «| αν «πο αρ «ο «ο απο «ο «ο «πο «πὸ s «ие s“ s «шә «шә 
«әб сш OF E E m E ЧЫР G G ЧЫР ФЕР ЧЫР «шә G «шә «шә G G G «шә «шә «жә «жә “ы: «ж» «нә ЧЫР ЧЫР < G G ЧЫР «и» «ив αν «| «ο αρ «αν s“ αν «| «ο απο «πὸ s s < s. G «но» “G G G G «и» «и» «не «и «нә «шә «шә 


SUBROUTINE READLCLLINES) 
INTEGER COUNT, I, NIX 
CHARACTERX4 BBLANK/'* '/, LLINESC8) 
DO 10 I-1,8 
LLINESCI) = BBLANK 
CONTINUE 
COUNT =0 
COUNT =COUNT +1 
IFCCOUNT.LT.3) GO TO 30 
WRITEC6,70) 
GO TO 50 
CONTINUE 
REWIND 5 
READC5,80, END=40, ERR=30)(LLINES(CJ),J=1,9) 
RETURN 
REWIND 5 
WRITEC6,60) 
GO TO 20 
CONTINUE 
STOP 


FORMATCIX,' WARNING: NULL STRINGS ARE NOT ALLOWED, ENTER CHARACTER 


X VALUES. ') 


FORMAT (///,5X,' PROGRAM TERMINATION - TWO NULL STRINGS ENTERED!') 


ҒОКМАТ(9А4) 
END 
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ОО ООО ООО 


еы со «тв «ше әже «и «ие «к» что «и что «ж. «тр «че чт s «е «то «и что «тр «тр “ше алы «лю» «то c «ж. s шз чыз аж» ты тыз шз сз тыз ο ο ο ο ошо еш» еше сы» «ә. «ж» «ж» сш» «ж» «жш» ош» 
ж“. «ж. ск» «же «к: ο ew ο ee «и» «то ех» a «и «чт «тр «и ee a «ят ew «тв что «те «тв чт «те «ят «те «те «но a ew «чт a a «и «те чт что s s s “ «и что «но «то «и «жо es s «ж» еже сш» «ню «же «жо оқ» өш» оқы оқ» «ж. πο at s s s 


SUBROUTINE META -- BY CALLING COMPRS AS A SUBROUTINE HERE, 


DONEPL HAS SUFFICIENT TIME TO FINISH; 


AND METAFILE NOT SAVED 


=o ab ab «ә» «тр «к: «к: «к: чи. «к. «к. «шо «ке «к» «кә «шә «к» «к. ei a «тр oii «к» «к. «к: πρ «| αὶ «αν αὶ αὶ αὶ αὶ αὶ «πὶ «πὶ, αὶ «αὶ «πὶ αὶ αὶ αὶ ab mb ab ab ab ab ab ab ab amb ab =b =b mb m 
- ο «тр «те «те «к: «те ab ab «к: «к» «к: чыз «ж» «ше ab s“. «тр «шә еш» «те «тр «#0 «тр «тр AD s. s. тыз =з тыз ЧЫР «тр «то «πὶ «ив «ие «и «πῶ «то αὶ «πὶ «πὶ, s. ab ab «и «тр «то ab ab mb ab чт «αὶ «πῶ «πὶ еш» «πὶ «πὰ αὶ өші cue ce c «ші» сш» еше 


10 


SUBROUTINE META 
DO 10 I-1,900000 
CONTINUE 

CALL COMPRS 
RETURN 

END 
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OTHERWISE COMPRS IGNORED 
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